Estimating the Instantaneous Frequency of Linear and Nonlinear Frequency Modulated Radar Signals—A Comparative Study

Automatic modulation recognition plays a vital role in electronic warfare. Modern electronic intelligence and electronic support measures systems are able to automatically distinguish the modulation type of an intercepted radar signal by means of real-time intra-pulse analysis. This extra information can facilitate deinterleaving process as well as be utilized in early warning systems or give better insight into the performance of hostile radars. Existing modulation recognition algorithms usually extract signal features from one of the rudimentary waveform characteristics, namely instantaneous frequency (IF). Currently, there are a small number of studies concerning IF estimation methods, specifically for radar signals, whereas estimator accuracy may adversely affect the performance of the whole classification process. In this paper, five popular methods of evaluating the IF–law of frequency modulated radar signals are compared. The considered algorithms incorporate the two most prevalent estimation techniques, i.e., phase finite differences and time-frequency representations. The novel approach based on the generalized quasi-maximum likelihood (QML) method is also proposed. The results of simulation experiments show that the proposed QML estimator is significantly more accurate than the other considered techniques. Furthermore, for the first time in the publicly available literature, multipath influence on IF estimates has been investigated.


Introduction
Recent years have seen dynamic development in radar instruments, including advancement in multiple-input multiple-output (MIMO) techniques, the pursuit to put quantum sensors into practice and the continual exploitation of even more ingenious modulation types. The last fact could not be overlooked from the context of electronic warfare (EW). This activity is aimed at gathering information about adversary radars, which is usually accomplished by highly sensitive interception receivers, specialized in collecting signals emitted by far, non-cooperative emitters. During peacetime, intercepted radar signals may be analyzed to find out the purpose and capabilities of the specific device. By contrast, during a conflict, the detection, identification and geolocation of enemy radar becomes far more important. The acquired intelligence can then be utilized to target or penetrate hostile air defense. Another relevant function of interception systems is supplying emitter data to databases used by other EW systems, such as radar warning receivers (RWRs).
Traditionally, electronic interception assets are divided into two distinct categories: electronic support measures (ESMs) and electronic intelligence (ELINT), where the former 1.
It is the first extensive comparison of five different instantaneous frequency estimation methods which are specific to radar signals recognition, which also accounts for various frequency modulation laws.

2.
We introduce an alternative IF estimation algorithm based on the generalized quasimaximum likelihood (QML) method with a novel signal model. 3.
Multipath propagation effects inherent to ESM analysis has been discussed and their influence on respective estimation methods is examined.
The paper is organized as follows. Section 2 presents the recognition system model as well as the assumed signal model. A multipath phenomenon description is introduced in Section 3. Section 4 gives the considered IF estimation methods. Section 5 shows the simulation experiments' results. Finally, Section 6 discusses the paper and presents the conclusions.

System and Waveforms
A typical automatic modulation recognition algorithm in an ESM system may be divided into the feature extraction stage and classifier stage, as shown in Figure 1. Here, we assumed a common situation, i.e., that one of the analyzed signal features is the IF. First, the signal was detected and then carrier frequency was estimated and removed. Hence, a complex envelope was determined. The time of arrival and pulse duration were assumed to be known. The AMC algorithm itself is not considered in this paper. different methods practiced for IF estimation with regard to linear and nonlinear frequency modulated radar signals.
The contributions of this paper are as follows: 1. It is the first extensive comparison of five different instantaneous frequency estimation methods which are specific to radar signals recognition, which also accounts for various frequency modulation laws. 2. We introduce an alternative IF estimation algorithm based on the generalized quasi-maximum likelihood (QML) method with a novel signal model. 3. Multipath propagation effects inherent to ESM analysis has been discussed and their influence on respective estimation methods is examined.
The paper is organized as follows. Section 2 presents the recognition system model as well as the assumed signal model. A multipath phenomenon description is introduced in Section 3. Section 4 gives the considered IF estimation methods. Section 5 shows the simulation experiments' results. Finally, Section 6 discusses the paper and presents the conclusions.

System and Waveforms
A typical automatic modulation recognition algorithm in an ESM system may be divided into the feature extraction stage and classifier stage, as shown in Figure 1. Here, we assumed a common situation, i.e., that one of the analyzed signal features is the IF. First, the signal was detected and then carrier frequency was estimated and removed. Hence, a complex envelope was determined. The time of arrival and pulse duration were assumed to be known. The AMC algorithm itself is not considered in this paper. A feature extractor is supplied with a discrete-time radar pulse corrupted with noise ( ) y n and sampled with rate s f . The signal model is then given by x n n a n e n where ( ) x n denotes the complex radar signal, n is the sample index, s T stands for sampling interval ( ) ( ) = s x n x nT , ε( ) n represents complex additive white Gaussian noise (AWGN) with variance 2 ε σ and ( ) a n with ( ) n φ are instantaneous amplitude and instantaneous phase, respectively. In reality, the pulse amplitude is not constant during the time, as both the rising and falling edge have some finite time. This phenomenon is usually caused by the transmitter, but the envelope can be also tapered deliberately in order to limit out-of-band emissions or suppress the sidelobes in matched filter. In the literature, the instantaneous amplitude of a pulse is frequently modeled by weighting functions such as the Tukey or raised cosine window. Here, similarly to [23], we adopted the Tukey window with the shape parameter α 0.05 T = . A feature extractor is supplied with a discrete-time radar pulse corrupted with noise y(n) and sampled with rate f s . The signal model is then given by where x(n) denotes the complex radar signal, n is the sample index, T s stands for sampling interval x(n) = x(nT s ), ε(n) represents complex additive white Gaussian noise (AWGN) with variance σ 2 ε and a(n) with φ(n) are instantaneous amplitude and instantaneous phase, respectively.
In reality, the pulse amplitude is not constant during the time, as both the rising and falling edge have some finite time. This phenomenon is usually caused by the transmitter, but the envelope can be also tapered deliberately in order to limit out-of-band emissions or suppress the sidelobes in matched filter. In the literature, the instantaneous amplitude of a pulse is frequently modeled by weighting functions such as the Tukey or raised cosine window. Here, similarly to [23], we adopted the Tukey window with the shape parameter α T = 0.05. This paper is focused on nonlinear frequency modulated radar signals. Along with LFM, it is one of the first conceived radar waveforms. NLFM's advantage over LFM comes from its naturally low sidelobes level, attainable without a need of mismatched processing, which is not the case for LFM. Thanks to this, the entire transmitter energy can be embraced. This is reflected in growing interest in NLFM, especially for synthetic aperture radars [24] or long-range surveillance radars [25]. The last decade brought rapid advancement in NLFM synthesis methods, with algorithms based on Bézier curves [26] or genetic optimization techniques [27]. The vast majority of already devised methods have one thing in common, namely that yielded IF functions are monotonic, symmetric with respect to time point corresponding to center of the pulse and the chirp rate is abruptly increasing towards both edges of the pulse. Different, i.e., nonsymmetrical IF generally would result in a deteriorated sidelobes level and Doppler tolerance [26] so it is hard to find such forms of NLFM in real-life applications. If we designate the constant-time variable as t, pulse duration as T, frequency deviation as ∆F, and IF as f (t) and further assume that pulse is defined in time period [−T/2, T/2], it may be concluded that the IF of the NLFM signal meets the following assumptions: The following NLFM signal synthesis methods were used in the remainder of this study: 1.
Based on the stationary phase principle (SPP). It is a primary design approach, originally devised by Fowle in the 1960s [28]. The IF is determined using its approximative relationship with inversed group delay, which in turn is obtained from the assumed prototype of energy spectral density. An interested reader can find a concise explanation of the method in [29] (pp. 87-89). Here, during the synthesis process, we picked Blackman-Harris due to its ability to achieve low sidelobes while preserving adequate range resolution [30].

2.
A design algorithm which describes the instantaneous phase of the NLFM waveform as a combination of the catenary and kappa curve [31]. In experiments, we used the optimization parameter α = 6 as defined in [31] and assumed that the maximum Doppler shift f D = 6 kHz.

3.
The NLFM waveform originally proposed for active sonar [32]. For our study, we skipped additional amplitude tapering and applied parameters of frequency modulation function α = 0.52 and γ = 1.47 as defined in [32]. Furthermore, the LFM was also employed in research for comparison. One can note that frequency variations in the NLFM waveform can have different sizes, depending on specific design requirements. For this reason, it is helpful to introduce the notion of nonlinearity factor, given by where f LFM (t) denotes the IF of the LFM signal with the same pulse duration, frequency deviation and modulation direction (here we assumed that both signals are upchirp) as the analyzed NLFM signal. Multiplication by eight is used to fit the coefficient values to the interval [0, 1], where 0 represent signals with purely linear frequency modulation and 1 describes signals which exhibit sharp frequency variations on the pulse borders (the waveform proposed by Price [29] (p. 92) may be used as an example of such). The nonlinearity factor can be simply seen as a value proportional to the area between the normalized IF of the NLFM and LFM signals, as illustrated in Figure 2. The comparison of different NLFM waveforms in the function of factors such as a time-bandwidth product (TBP) is further complicated by the fact that for the NLFM signal, the frequency deviation is not equal to the effective signal bandwidth B (generally, deviation is larger). For this reason, in experiments we assumed that the instantaneous frequency took the values from the interval f (t) ∈ [−∆F/2, ∆F/2], which are adjusted in order to fit the required bandwidth. Meanwhile, the bandwidth is interpreted here as the frequency span which contains the majority of the signal energy. To facilitate this, the in-bound energy ratio parameter (IBER) may be used, introduced in paper [33] and defined as where |X( f )| represents the amplitude spectrum of a signal. IBER = 0.9 is used in simulation experiments. Finally, all radar waveforms employed in the study are summarized in Table 1.   Side effects which come from the multipath were thoroughly investigated in the context of global navigation satellite systems (GNSS) and wireless communications. Nonetheless, multipath is averse to the ESM systems as well, while its most recognizable outcome is the deterioration of direction-finding algorithms. Until recently, only a few studies have taken up this subject, such as [34] (p. 266). Meanwhile, significant contribution was made by Robertson in her book [1]. Here, we adopt a simplified description of the phenomenon, yet sufficient to evaluate the considered algorithms. Inspired by the profound DMCM model in [35], if we assume that each signal cluster consist of one propagation path, both the radar and ESM platform are stationary and no local scattering is present, the pulse affected by the multipath y M (t) can be represented as a vector sum of the direct-path pulse y D (t) and its Z reflected copies y R (t): where λ < 1 denotes the reflection coefficient which is related to the attenuation level at the reflection surface and depends on the grazing angle, the polarization of the reflected signal and the admittance of the reflecting medium [1]. Further, ∆t is the time lag between the direct-path signal and reflected signal (for the intra-pulse distortions to occur ∆t < T). Clearly, reflected rays will be delayed with respect to the direct-path signal, making a detour of ∆R =c ∆t, where c stands for speed of light in vacuum. Next, ζ denotes the amplitude attenuation factor accounting for path loss, whereas power attenuation is given by P dB = 20 logζ [36]. Following power delay profiles remarked in [35], here we assume that power decreases linearly with respect to distance [36] and P dB = −0.03∆R. Next, ∆φ represents the phase difference between the direct and reflected signal, and z is a number of the respective reflected pulse.
The multipath phenomenon is commonly described by statistical means. In simulation studies, we assumed that the phase has a uniform distribution as propagation delays are usually much greater than the carrier signal period, so all phase values are equally likely. Next, for a given wavelength, the reflection coefficient strongly depends on the reflection geometry and scattering medium [1]. As we do not model any specific reflectors nor spatial scenario, we adopt that λ is governed by the uniform distribution. Lastly, we assume that time delays ∆t follow Gamma distribution, as found in [36]. This observation relates to GNSS, albeit in the public domain there are not adequate studies which specifically concern ESM. A diagram of the multipath phenomenon for a single reflection is depicted in Figure 4a. It should be noted that the phase difference between consecutive signals determines whether interferences will take a destructive or constructive form, as can be inferred from Figure 4b.
tors nor spatial scenario, we adopt that λ is governed by the uniform distribution. Lastly, we assume that time delays t Δ follow Gamma distribution, as found in [36]. This observation relates to GNSS, albeit in the public domain there are not adequate studies which specifically concern ESM. A diagram of the multipath phenomenon for a single reflection is depicted in Figure  4a. It should be noted that the phase difference between consecutive signals determines whether interferences will take a destructive or constructive form, as can be inferred from Figure 4b.

Considered IF Estimation Methods
In radar AMC algorithms, there are two widely used IF estimation techniques: based on phase finite differences and time-frequency distributions. Here, we scrutinize both approaches, considering five distinct estimation methods. Three of them are variants of phase differentiation techniques (Section 4.1), two employ TFDs (Sections 4.2 and 4.3) and the remaining is based on the generalized QML method (Section 4.4).

Phase Finite Differences
The IF of a signal is a derivative of its instantaneous phase ( ) t φ with respect to time, i.e., whereby for a discrete-time analytic, a signal's instantaneous phase can be easily determined as an inverse tangent of the ratio of quadrature and the in-phase component.
Phase function should be unwrapped in prior, i.e., its absolute jumps towards π should be replaced by their 2π ± complements. After the discrete, unwrapped instantaneous phase of the intercepted signal ( ) n φ is obtained, the IF estimate  ( ) f n can be easily calculated using a popular numerical method for solving differential equations, namely by a finite difference of a form

Considered IF Estimation Methods
In radar AMC algorithms, there are two widely used IF estimation techniques: based on phase finite differences and time-frequency distributions. Here, we scrutinize both approaches, considering five distinct estimation methods. Three of them are variants of phase differentiation techniques (Section 4.1), two employ TFDs (Sections 4.2 and 4.3) and the remaining is based on the generalized QML method (Section 4.4).

Phase Finite Differences
The IF of a signal is a derivative of its instantaneous phase φ(t) with respect to time, i.e., whereby for a discrete-time analytic, a signal's instantaneous phase can be easily determined as an inverse tangent of the ratio of quadrature and the in-phase component. Phase function should be unwrapped in prior, i.e., its absolute jumps towards π should be replaced by their ±2π complements. After the discrete, unwrapped instantaneous phase of the intercepted signal φ(n) is obtained, the IF estimatef (n) can be easily calculated using a popular numerical method for solving differential equations, namely by a finite difference of a form The aforementioned method is termed as backward finite difference (BFD). It was referred to both in the ESM literature [34] (p. 270) and afterwards exploited in numerous AMC algorithms [7,8,16,37] for feature extraction of a signal.
Similarly, after applying a minor amendment in the above formula, a central finite difference (CFD) estimator is yielded. It is unbiased, does not introduce group delay for LFM signals [18] and can be noted aŝ Both mentioned IF estimators exhibit high variance while dealing with distinct noise. There are a couple of solutions to limit the variance. One of them is to average or "smooth" successive phase differences in the moving window. The estimator proposed by Kay (hereinafter referred to as KAY) may serve here as an example [38] with an averaging window of a length N − 1, formed by Due to windowing, the variance of the estimate is reduced by approximately N/6 [18]. In simulation experiments, we adopted N = 16, which assures a substantial error reduction. The application of a considerably longer window is inadvisable since the estimate will be degraded if rapid frequency variations are present within the window [18].

Data Driven Pseudo-Wigner Distribution
Another prevalent method of IF estimation is to extract the ridge from the timefrequency image of a signal [18,19]. Authors of the AMC algorithms also draw from this approach, using different forms of Wigner-Ville distribution (WVD), as can be seen in [11,39]. Next, the IF estimator in consideration belongs to this category. In the first step, a time-frequency image is obtained through pseudo-Wigner distribution (PWD) with an adaptive window length. The discrete PWD is described by the formula [40] where N is a positive integer, k is the frequency index and w(m) denotes the real-valued window of length 2N. The idea behind the algorithm is to select length of the PWD window independently for each time instant, in order to reach best bias-variance tradeoff. An interested reader may find its details in [40]. In simulations, we used a rectangular window with lengths ranging from 10 to 128 samples. Wigner distribution has the best energy concentration among quadratic TFDs [41] (p. 585). PWD also exhibits good concentration and therefore the algorithm assumes that for a given time instant, the IF estimate can be simply obtained from the peak of the time-frequency distribution, i.e.,f where Q k denotes set of frequency samples obtained for an optimum window length. The same IF estimation strategy was employed in a notable AMC algorithm developed by Lunden [7] and later by the authors of [42]. Similar to [7], we adopted additional median filtering of the IF estimate with a window size of 10 samples. This removes "spikes" which emerge in the estimate for low signal-to-noise ratios (SNRs) [43].

Choi-Williams Distribution Image
Another TFD often used in radar AMC algorithms [7][8][9] is Choi-Williams distribution (CWD), also known as exponential distribution. Its continuous form is given by [44] (p. 92) where τ denotes the lag variable and scaling factor σ controls the tradeoff between crossterms suppression and frequency resolution. For experiments, we picked a value σ = 0.05, likewise to [7]. Here, we determine CWD by multiplying its kernel g(ν, τ) = exp(−ν 2 τ 2 /σ) with WVD result in the ambiguity domain [44] (p. 245), which may be represented as where ν denotes Doppler shift, A x (ν, τ) stands for symmetric ambiguity function and W WV (t, f ) represents the WVD. CWD's inhering ability to attenuate the cross-terms is at the expense of worsened energy concentration. For this reason, in the case of a CWD-based estimator, we do not exploit the time-frequency ridge, but the morphological IF extractor is in its place. Various AMC algorithms [9,39] have used mathematical morphology before. Here, we propose an image processing algorithm based on one described in [8]. The absolute value of the CWD image forms an input image CWD(i, l), where i, l denotes pixel indices. Then, adaptive binarization is performed as follows: a.
Calculate the normalized grayscale image with the values from the interval [0, 1], as Set values of all pixels which meet the condition G(i, l) ≤ V to 0; d.
Calculate the Shannon entropy of the resultant image where probability associated with the given pixel value g r is estimated using the 256 bins histogram and R denotes the total number of gray levels; e.
Decrease the threshold value by specified step V = V − δV; in this paper δV = 0.005 is used; f.
compute final binary image B(i, l), given by After binarization, some isolated noise remains in the image. In order to limit the noise and emphasize the main trend in the IF, the image is additionally treated with morphological opening followed by closing, with a structuring element of increasing size. We employed a square with a size 3 · 3-5 · 5. In the end, the pretreated image may be interpreted again in the time-frequency domain, by simply associating individual pixels with the corresponding time and frequency values B(n, k).
Lastly, an instantaneous frequency estimate is calculated for each time instant as a median of all K frequency values corresponding to non-zero pixelŝ f (n) = median[B(n, 0), B(n, 1), . . . , B(n, K)].
The workflow of the CWD-based estimator is depicted in Figure 5. Two subplots on the left-hand side and W CW (·) are related to the transforms described by formulae  (14) and (15). The three remaining subplots illustrate the operation of the morphological ridge extractor. pixels with the corresponding time and frequency values ( , ) B n k .
Lastly, an instantaneous frequency estimate is calculated for each time instant as a median of all K frequency values corresponding to non-zero pixels The workflow of the CWD-based estimator is depicted in Figure 5. Two subplots on the left-hand side and ( ) CW W  are related to the transforms described by formulae (14) and (15). The three remaining subplots illustrate the operation of the morphological ridge extractor.

Generalized QML Method
Several authors have remarked that the NLFM waveform may be accurately modeled as a polynomial-phase signal (PPS) [20,24,33]. Thus, in this work we adopt the quasi-maximum likelihood (QML) IF estimator, originally proposed by Djurović and Stanković [45] for PPS signals. However, the initial experiments showed that fast frequency variations on pulse borders prevent signals from being accurately modeled by polynomials only. Consequently, in order to improve the convergence of the model at the borders, we enriched it with tangent component. Finally, a generalized IF model is given by

Generalized QML Method
Several authors have remarked that the NLFM waveform may be accurately modeled as a polynomial-phase signal (PPS) [20,24,33]. Thus, in this work we adopt the quasi-maximum likelihood (QML) IF estimator, originally proposed by Djurović and Stanković [45] for PPS signals. However, the initial experiments showed that fast frequency variations on pulse borders prevent signals from being accurately modeled by polynomials only. Consequently, in order to improve the convergence of the model at the borders, we enriched it with tangent component. Finally, a generalized IF model is given by where a i denotes the respective coefficients of the M-th order polynomial. In experiments, we adopted M = 4 in order to avoid overfitting of the model which might arise for higher order polynomials. Furthermore, it is assumed that a 1 = a 3 = 0 to account the fact that frequency modulation law (2) is symmetric. The concept of this method is to compute a series of IF estimates using discrete shorttime Fourier transform (STFT) W SF(h) (t, f ), determined with different window lengths from set H. The estimate which maximizes the QML function is considered to be the most accurate. In research, we applied a set of rectangular windows of widths h i = 2 (i+1) samples, where i means integer from the interval [1, I], and maximum window length is limited by the number of available signal samples h I ≤ N.
The algorithm is summarized below, whereby each step is performed for each window length h ∈ H (for brevity we use continuous time and frequency arguments): 1.
Short-time Fourier transform W SF(h) (t, f ) is obtained.

2.
The initial IF estimate f (t) is obtained from the ridge of the spectrogram 3.
Coefficients of the signal model (17) are estimated using least-squares polynomial regression function F(·) F(â 0 ,â 1 ,â 2 ) = argmin 4. Particular frequency, and then phase is reconstructed on the base of estimated coefficientsâ i , asf Finally, among all calculated estimatesf h (t) we select those that is maximizing the QML function of a form [45] Moreover, in simulations, we adopted a slight amendment to the algorithm. It stems from the fact that absolute values of the IF estimate purse to the infinity on waveform edges. Hence, values of 0, 2%T f s boundary samples (where · denotes ceiling function) on the left and right edge of the pulse are substituted by values of succeeding and proceeding samples adequately.

Results and Discussion
This section presents the results of conducted experiments. The necessary software was written and executed in MATLAB 2018 and code is available to the reader on request. Algorithms were evaluated using set of complex, baseband test signals buried in noise. The sampling frequency was set to f s = 100 MHz. The number of available signal samples was determined by the signal length. Pulse duration and bandwidths were random variables to reflect the fact that signals intercepted by the intelligence receiver have varying parameters. The distribution of parameters applied in simulation trials is given in Table 2. We consider SNR values from a range typical to ESM systems, i.e., −5 ÷ 30 dB, whereby the SNR is defined as 10 log(P x /σ 2 ε ) and P x means the average power of a signal estimated using the mean square value. Symbol U(a, b) signifies that the corresponding parameter is evenly distributed within range a to b, Γ(ς, θ) means that it takes a value from gamma distribution with shape ς and scale θ and [a, b] denotes any integer value within given interval. The accuracy of the algorithms is established using the mean squared error (MSE) evaluated on 250 Monte Carlo trials, i.e., where N trial denotes the number of Monte Carlo runs,f t (n) is the IF estimate with a length of N samples obtained in t-th trial, and f (n) stands for the true IF law of one of the evaluation signals introduced in Section 2.

General Accuracy
In first stage, we compared the accuracy of respective methods, when dealing with LFM as well as various forms of NLFM signals. Figure 6 presents the results obtained for the NLFM A waveform. Our generalized QML algorithm proved to be the most precise among compared methods, across the entire range of SNR, even for its low values (it stems from the fact that STFT is relatively robust to the noise influence despite being biased). The accuracy of both the CWD-based method and data-driven PWD are comparable, albeit the former has a disadvantage in being also reliant on the image processing stage performance. It explains the increase in CWD error, observed for higher SNR values.
For clarity, the results obtained for different waveforms are grouped in Table A1 in Appendix A. However, they follow the same trend. Exemplary instantaneous frequency estimates obtained using respective algorithms for SNR = 5 dB are shown in Figure 7. Exemplary instantaneous frequency estimates obtained using respective algorithms for SNR = 5 dB are shown in Figure 7.

Accuracy for Different Nonlinearity Levels
When examining Figure 7, one can discern that for the NLFM signal, the significant part of the error results from the poor accuracy of the estimate on pulse borders. This is reaffirmed in next experiment, whose results are shown in Figure 8a. In this case, MSE was determined entirely on pulse edges. It turns out that the error value is increased by approximately 1 dB, regardless of the applied IF estimation method. In other words, although all estimators perform quite well when the frequency is shifting linearly, their resolution turns out to be insufficient when varying chirp rate come into play (notwithstanding that both PWD and QML employ windows of adaptive size). This is illustrated in Figure 8b. In practice, it may result in the misclassification of the NLFM as an LFM in the AMC algorithm. Exemplary instantaneous frequency estimates obtained using respective algorithms for SNR = 5 dB are shown in Figure 7.

Accuracy for Different Nonlinearity Levels
When examining Figure 7, one can discern that for the NLFM signal, the significant part of the error results from the poor accuracy of the estimate on pulse borders. This is reaffirmed in next experiment, whose results are shown in Figure 8a. In this case, MSE was determined entirely on pulse edges. It turns out that the error value is increased by approximately 1 dB, regardless of the applied IF estimation method. In other words, although all estimators perform quite well when the frequency is shifting linearly, their resolution turns out to be insufficient when varying chirp rate come into play (notwithstanding that both PWD and QML employ windows of adaptive size). This is illustrated in Figure 8b. In practice, it may result in the misclassification of the NLFM as an LFM in the AMC algorithm.
(a) (b) Figure 8. (a) The mean squared error achieved by considered IF estimation methods, evaluated on pulse edges only. The error was determined for time samples corresponding to 7% of pulse duration for both pulse ends. The NLFM A waveform was concerned, as on Figure 6; (b) illustration of how an inaccurate IF estimate may cause misclassification in the recognition system.
To better understand this effect, we performed additional calculations of the estimation error in a function of the nonlinearity factor. The results are presented in Figure 9. The error was determined for time samples corresponding to 7% of pulse duration for both pulse ends. The NLFM A waveform was concerned, as on Figure 6; (b) illustration of how an inaccurate IF estimate may cause misclassification in the recognition system.
To better understand this effect, we performed additional calculations of the estimation error in a function of the nonlinearity factor. The results are presented in Figure 9. ure 8. (a) The mean squared error achieved by considered IF estimation methods, evaluated on pulse edges only. The or was determined for time samples corresponding to 7% of pulse duration for both pulse ends. The NLFM A wavem was concerned, as on Figure 6; (b) illustration of how an inaccurate IF estimate may cause misclassification in the ognition system.
To better understand this effect, we performed additional calculations of the estimation error in a function of the nonlinearity factor. The results are presented in Figure 9. The noteworthy property of finite difference estimators is their invariance to modulation form, despite being far less accurate when compared with other methods. The same relates to the PWD and CWD estimator. Concededly, both techniques cope slightly better with modulation laws which are close to linear (it emerges that both TFDs exhibit The noteworthy property of finite difference estimators is their invariance to modulation form, despite being far less accurate when compared with other methods. The same relates to the PWD and CWD estimator. Concededly, both techniques cope slightly better with modulation laws which are close to linear (it emerges that both TFDs exhibit better energy concentration for the LFM), but the decrease in their accuracy for NLFM is minimal.
The proposed parametric QML estimator again revealed itself to be the most precise. The devised frequency modulation model (17) nicely fits the IF-law of NLFM waveforms for nonlinearity factors smaller than 0.4. Above this value, the estimation error rises considerably, but still remains lower in comparison to other methods (it is to be noted that NLFM waveforms of high η will have a strongly diffused spectrum which refrains radar designers from implementing such waveforms in real-life). Still, it is covetable to find another regression model which would better fit the modulation laws of applicable NLFM waveforms with an η value equaling at least 0.6.

Accuracy in Multipath Conditions
Available reports state that the multipath in ESM receivers appears primarily by quasirandom fluctuations in the pulse amplitude [1]. At the same time, the IF is traditionally considered as a reliable signal feature, which may even carry unintentional modulation on pulse (UMOP) which is characteristic to a certain emitter and thus usable for fingerprinting [46]. Hence, in the next step we surveyed the multipath impact on respective IF estimators. In simulations, we represented the multipath signal as a sum of the discretetime complex radar signal and its delayed replicas with accordance to formula (5), whereas respective parameters of the model were considered to be random variables with distributions given in Table 2. Complex AWGN noise samples were added to the composite signal, not to respective pulses. Figure 10 illustrates IF estimates obtained using different methods, in the case of a single reflection with an amplitude amounting 0.9 of that of the direct path component is present in the signal. Such a high value of the alternative path component may be assumed as a strong reflection, and may occur for small grazing angles. the discrete-time complex radar signal and its delayed replicas with accordance to formula (5), whereas respective parameters of the model were considered to be random variables with distributions given in Table 2. Complex AWGN noise samples were added to the composite signal, not to respective pulses. Figure 10 illustrates IF estimates obtained using different methods, in the case of a single reflection with an amplitude amounting 0,9 of that of the direct path component is present in the signal. Such a high value of the alternative path component may be assumed as a strong reflection, and may occur for small grazing angles. It can be noted that for a considered time, the phase and amplitude relationship between the direct-path and reflected pulse, the envelope of the resultant pulse (seen by It can be noted that for a considered time, the phase and amplitude relationship between the direct-path and reflected pulse, the envelope of the resultant pulse (seen by the ESM system) is severely distorted by destructive interference. For time instants where the amplitude almost completely fades, the estimation error clearly has the biggest value. As it is distinct from other methods, QML does not produce error peaks as it is the only parametric estimator in consideration. To better examine the outcomes of multipath, we performed two simulations aiming at assessing the increase in estimation error for the multipath scenario, compared to a signal free of reflections. The IF of the direct-path pulse was used as the reference for MSE calculation (in fact it is imprecise while true IF is also altered by phase shifts related to the multipath). The results are presented in Figure 11. As it is distinct from other methods, QML does not produce error peaks as it is the only parametric estimator in consideration. To better examine the outcomes of multipath, we performed two simulations aiming at assessing the increase in estimation error for the multipath scenario, compared to a signal free of reflections. The IF of the direct-path pulse was used as the reference for MSE calculation (in fact it is imprecise while true IF is also altered by phase shifts related to the multipath). The results are presented in Figure  11.
(a) (b) Figure 11. Increase in estimation error in multipath conditions. In order make the results independent of different performance of the estimators at the pulse edges, only the central part of the pulses was taken into account, i.e., 7% of samples on both pulse ends were excluded from the calculation of the MSE. The NLFM A waveform was used. (a) Error increase as a function of SNR. Simulation parameters from Table 2 are applied. (b) Error increase for a single reflection with fixed ratio of direct pulse to reflected pulse and determined SNR = 15 dB.
Perhaps unsurprisingly, the estimation error grows for all methods, albeit to different extents. Referring to Figure 11a, the error growth is rather unsubstantial for finite phase difference methods, and even declines when the SNR falls below 10 dB. This is due to fact that those estimates are already noise-like, thus MSE approaches its maximum. One can note that KAY is noticeably more vulnerable to the multipath in comparison to BFD and CFD. Although all three estimators rely on the signal envelope, the former Figure 11. Increase in estimation error in multipath conditions. In order make the results independent of different performance of the estimators at the pulse edges, only the central part of the pulses was taken into account, i.e., 7% of samples on both pulse ends were excluded from the calculation of the MSE. The NLFM A waveform was used. (a) Error increase as a function of SNR. Simulation parameters from Table 2 are applied. (b) Error increase for a single reflection with fixed ratio of direct pulse to reflected pulse and determined SNR = 15 dB. Perhaps unsurprisingly, the estimation error grows for all methods, albeit to different extents. Referring to Figure 11a, the error growth is rather unsubstantial for finite phase difference methods, and even declines when the SNR falls below 10 dB. This is due to fact that those estimates are already noise-like, thus MSE approaches its maximum. One can note that KAY is noticeably more vulnerable to the multipath in comparison to BFD and CFD. Although all three estimators rely on the signal envelope, the former "smooths" its values, which results in error amplification when the amplitude exhibits a longer droop.
In the case of TFD-based methods, degradation in accuracy is more evident, which can be noted both from Figure 11a,b. This arises from the cross-terms originating from reflected pulses, as previously remarked in [20]. However, CWD and QML perform slightly better as-when compared to PWD-they have a better ability to suppress the cross-terms, which are further attenuated by the morphological IF extractor and regression procedure, respectively. Nevertheless, dropping performance of the QML is noticeable as the method is designed for slow variations in the instantaneous amplitude. It remains for further study of the improvement of the QML that can be used with high accuracy for multipath signals.

Conclusions
The choice of proper IF estimator for radar AMC algorithm is application-specific. If a certain ESM system copes with signals with high to moderate SNR values (above 10-15 dB), finite difference-based methods are sufficient. This approach clearly benefits in simplicity-calculations may be performed in a moving window, which facilitates hardware implementation. Among three analyzed estimators of that class, i.e., BFD, CFD and KAY, the latter achieves the lowest variance, which for high SNR approaches that characteristic to TFD-based methods. Still, if the intercept system needs to deal with weaker signals and extra processing power is available, TFD-based algorithms can be the solution, since they are relatively robust to the noise influence.
The highest accuracy can be obtained using parametric estimators. However, they require a priori knowledge about the modulation type present in the signal. This results in three-step signal analysis: firstly, the initial IF estimate may be obtained using general estimators (e.g., CFD), then the AMC algorithm (equivalently human operator) senses the modulation type, and finally the target estimation method suitable for a certain waveform may be employed to refine the initial estimate. Hence, parametric methods are particularly useful when real-time operation is not required, the situation characteristic, e.g., to classical ELINT systems.
The proposed generalized QML method is able adequately reconstruct nonlinear frequency modulation embedded in the waveform. It can be further improved by developing another IF model, e.g., the tansec one described in [47]. There are couple of other prominent estimation methods, such extended generalized chirp transform (EGCT). Still, if one has to deal with a real-time requirement specific to ESM, their applicability is limited due to high computational burden. In this case, a generalized QML algorithm can be used in place.
The provided multipath simulations show that respective IF estimates may be corrupted by reflections. This is accompanied by severe distortions in instantaneous amplitude. What should be regarded in AMC algorithms is that different estimators have different tolerance to this phenomenon. It is an interesting prospect to apply techniques for IF estimation for multicomponent signals [44] (pp. 588-596) in this case, as a pulse affected by the multipath may be considered as a multicomponent signal. All that notwithstanding, we are going make an effort to modify the QML so that it can work more accurately for multipath signals.
Although our simulations concern mainly the NLFM waveform, the presented results may also be referred to different kinds of frequency-modulated radar signals, such as polyphase-coded frequency modulation (PCFM). The introduced notion of the nonlinearity factor may be useful in further research on the synthesis and recognition of NLFM signals. Table A1. Comparison of IF estimators' accuracy, for various linear and nonlinear frequency modulated waveforms, in a function of SNR. It is worth noting that the accuracy of finite difference methods is actually independent of the signal form.