Automatic Classification of Frequency-Modulated Radar Waveforms Under Multipath Conditions

Techniques for automatic modulation classification (AMC) of radar signals are crucial for spectrum sensing, passive bistatic radars, and electronic intelligence (ELINT). Most of the existing AMC algorithms were evaluated using solely synthetic data. Meanwhile, signals intercepted in real-life scenarios are extra distorted due to the multipath propagation. This article presents the novel, pattern recognition AMC method for frequency-modulated radar signals, with improved resistance to noise and multipath. This has been achieved by a multistage feature selection process, involving over a dozen of popular radar signal metrics. Based on the feature selection results, the two advanced waveform features were utilized in the designed algorithm, i.e., quasi-maximum likelihood (QML) instantaneous frequency (IF) estimate and fractional Fourier transform (FrFT) profile. The proposed AMC method has been evaluated using both synthetic and real data. The obtained results show that proposed classification framework achieves an overall accuracy of 93.6% for a set of real-life signals acquisitions, corrupted both by noise and multipath influence.


I. INTRODUCTION
A LGORITHMS for distinguishing unknown, intrapulse modulation embedded in the intercepted radar signals have been intensively developed over the last two decades. Such techniques are known in literature as automatic modulation recognition/classification (AMR/AMC). They are useful in multiple civil and military applications, including spectrum monitoring, passive bistatic radars, electronic support (ES), electronic intelligence (ELINT), and radar warning receivers (RWRs). For example, in ELINT/ES systems, prior modulation recognition facilitates deinterleaving of pulse trains originating from different emitters and recognition of the specific type of radar among numerous existing types. During the pulse technical analysis, some of the radar capabilities may also be inferred on the grounds of applied waveform. The extensive research of AMC methods was initiated in early 2000s [1]. The breakthrough was provided with the work of Lundén and Koivunen [2], where the first comprehensive classification system, suitable for most widely used radar waveforms, was proposed. Subsequent studies were aimed to In fact, the vast majority of the existing AMC algorithms are pattern recognition machine learning methods with different complexities of feature extraction and classifier stages. These algorithms proved to be quite effective for noise-buried synthetic signals, yet they have several drawbacks in common. 1) They have not been validated using real-life signals.
2) Usually, it is assumed that the signal is corrupted by the additive white Gaussian noise (AWGN) only. The primary goal, consequently, is to design the algorithm, which will remain accurate for the lowest possible signal-to-noise ratio (SNR). Meanwhile, in classic ELINT scenarios, received signals are easily distinguishable from the noise. On the other hand, the multipath propagation-phenomenon that may introduce severe distortions in the intercepted signal-has been scarcely noticed so far. 3) Many complex signal processing methods, such as time-frequency distributions, have been already applied in AMC. However, ELINT algorithms should be fit for real-time execution, what is barely considered.
For instance, few hardware AMC implementations reported [16], [17], [18], had to utilize simple features extracted from instantaneous signal properties (ISPs) in order to limit computational burden. This article presents new AMC algorithm, which solves the aforementioned problems. Instead of customary application of some arbitrarily chosen features, we compared a large set of candidate features and examined them in the aspect of tolerance to noise and multipath distortions. This allowed to design straightforward, yet accurate classification algorithm.
It should be pointed that few other attempts to mitigate the multipath outcomes has been already made both in the field of radar and communications signal processing, e.g., channel equalization [15], [19] or direct signal extraction [20]. However, we assumed that thorough feature selection will be sufficient to avoid further preprocessing.
Our current research is limited to two radar modulations: presumably, the most common [21] linear frequency modulation (LFM) and more advanced nonlinear frequency modulation (NLFM) with a few of its existing variants. Our interest in NLFM modulation stems from its increasing popularity [22]. In contrast to plain chirp, NLFM waveforms do not require weighting and thus do not induce mismatched filtering. Furthermore, unlike phase shift keying (PSK), NLFM has not been regularly considered in AMC so far.
The remainder of this article is outlined as follows. First, a recognition system and utilized waveforms are introduced in Section II. Then, Section III explains the multipath propagation model used throughout research. Next, Section IV gives a set of ten signal features under consideration, including time, time-frequency, and spectral representations of a signal typical for radar AMC. Afterward, Section V describes a threestage feature selection process. Section VI presents the utilized naïve Bayes (NB) classifier with its decision rule. Next, the designed AMC algorithm is verified using simulation and reallife signals, and the results are given in Section VII. Finally, this article is concluded in Section VIII.

II. SYSTEM AND WAVEFORMS
In this section, we introduce the requirements for designed AMC algorithm, present its role within the recognition system, and define five radar waveforms considered within the research.

A. Recognition System
The typical waveform classification system comprises three major parts, as shown in Fig. 1. The front-end receiver intercepts the high-frequency radar signal x RF (t), where t is the continuous-time variable, and then estimates its carrier frequency. Next, during complex downconversion, carrier is removed and signal is sampled at f s = 1/T s . Thereafter, signal is preprocessed in order to isolate respective pulses, so that time of arrival (TOA) and pulsewidth (PW) of each pulse are estimated-this if often done using adaptive thresholding methods. Finally, we can assume that we deal with consecutive N samples of a complex signal x(n), which is a single radar pulse s(n), transmitted over AWGN channel, i.e., where n denotes the sample index, a(n) is the instantaneous amplitude, and φ(n) stands for instantaneous phase (IP). The noise component ε(n) has the variance of σ 2 ε and the SNR is defined as 10log(P s /σ 2 ε ), where P s represents the signal's average power.
The signal x(n) is fed into a feature extractor and then classifier, which jointly constitutes the AMC algorithm. This prerequisites one-time feature selection and classifier training, whereby both are done offline. At the output, the recognition system yields the information that x(n) embeds one of the known modulations in this research, NLFM or LFM.

B. Considered Waveforms
In previous AMC studies which recognized NLFM as a separate class, e.g., [16], [23], only a single nonlinearity model has been considered. Meanwhile, no specific form of this waveform dominates. For this reason, in order to achieve proper generalization by the designed algorithm, during the research, we applied four different NLFM models, derived from the most popular synthesis methods reviewed in [22]. These waveforms are gathered in Table I, where f (t) = φ ′ (t)/2π signifies instantaneous frequency (IF), η denotes the nonlinearity factor [24], and PSL stands for peak sidelobes level.
All waveforms above have f (t) given by monotonic and odd function. This complies with the majority of NLFM synthesis algorithms [24], [29], [30]. Already existing AMCs did not account for this aspect during NLFM modeling. This might result from a gap between radar and reconnaissance communities.
It should be noted that in case of the NLFM, the discrepancy between frequency deviation δF and bandwidth B may be substantial. To prevent notion ambiguities, we adopted energetic definition of bandwidth using in-bound energy ratio (IBER) defined in [31] and scaled IF appropriately [22] to reach desired IBER = 0.9. Furthermore, to account for finite rise and fall time of a pulse, x(n) was weighted by discrete Tukey window to shape both slopes to 250 ns length.

III. MULTIPATH MODEL
In practice, the academic assumption that AWGN is the one factor which impairs the reception generally does not hold true. It should also be accounted that signal is scattering from different objects along its path, e.g., buildings, foliage, or ground. In turn, the multicomponent signal induced in the antenna is the product of the desired direct path signal and unwanted, delayed rays as depicted in Fig. 2. The multipath distortions in question have been carefully examined in a field of wireless communications and global navigation satellite systems (GNSS). The corresponding studies for electronic warfare were limited until the monograph [32]. Recently, Kong et al. [14], [15] introduced the multipath fading in a context of radar AMCs, giving momentum to this problem.
In the following, we present our multipath channel model utilized throughout feature selection and next in simulation experiments. Let us recall that the compound signal is the incoherent sum of s(n) and its R delayed replicas. Then, (1) can be rewritten in continuous time domain as where ξ signifies reflection coefficient, ζ represents amplitude attenuation factor, δt denotes time lag, and δφ is phase offset, whereas all parameters are specific for individual reflection. The underlying problem is to estimate these quantities. In communications, it is common to describe fading channel with its impulse response or power delay profile (PDP), whereby PDP can be determined in field measurements. Analogical characteristics for radar are not publicly available. Hence, we adopted qualitative approach and predetermined certain distribution for each parameter based on [24], [32], and [33].
The reflection coefficient ξ describes attenuation level at the reflection surface and depends on the grazing angle, admittance of the reflecting medium, and polarization of the incident signal [32], [34]. Since no geometrical dependencies are modeled, it was assumed that the coefficient takes values from the uniform distribution ξ ∼ U (0, 1).
Next, ζ provides amplitude scaling resulting from path loss, which is related to power attenuation as P att = 20log(ζ) [33]. The free-space path loss (FSPL) may be estimated using the formula [35] FSPL dB = 32.4 + 20log 10 (f c ) + 20log 10 (D km ) where f c denotes the carrier frequency expressed in MHz and D km is the distance from the transmitter, given in kilometers. One can find that for a typical microwave band used in radar and a distance of several dozens kilometers, FSPL function is locally linear. Given that, we performed simple linear regression on experimental PDP from [36], obtaining the following power attenuation function: where ∆ = cδt denotes detour and c is the speed of light. It can be inferred that multipath inference arises when time delay between direct path and reflected pulse δt < T . This effect is mostly visible by deformation of the pulse envelope, as illustrated in Fig. 3.
Based on [33], it has been supposed that ∆ is governed by gamma distribution with shape k Γ = 2.6 and scale θ Γ = 129, as depicted in Fig. 4. The model assumes that the phase shift of the baseband signal x(t) is distributed according to δφ ∼ U (0, 2π), as propagation delay is usually many times higher comparing to the carrier signal period. Finally, the number of paths R constitutes simulation parameter.

IV. FEATURE EXTRACTION
In this section, a set of considered features are presented. During this study, we examined many more features, including higher order moments and cumulants, Rényi entropy, and TFR variance. Nonetheless, in order to keep the article concise, we report only those features that were found informative enough after the initial selection. The section covers five feature categories specific to radar AMC. First, the features based on instantaneous signal properties are discussed. Next, a single feature based on the moment of the complex envelope is given. Then, the features derived from time-frequency distributions are proposed. After that, the features based on power spectral density (PSD) of a signal are described. Finally, new feature based on fractional Fourier transform (FrFT) is introduced. The ith feature is denoted as γ i .

A. Instantaneous Signal Properties
The features derived from instantaneous amplitude, frequency, and phase are frequently applied for AMC due to ease of their implementation. They were made popular by Azzous and Nandi [37] for communication signals recognition. In this research, we utilized IP and IF extracted using common phase differentiation methods, but also more advanced quasi-maximum likelihood (QML) algorithm. Instantaneous amplitude was omitted, since it is indiscriminative for LFM/NLFM pair.
1) ISP Derived From the Complex Envelope: The IP of s(n) may be estimated as where i(n) and q(n) denotes the in-phase and quadrature components of the recognized signal x(n), respectively,• signifies the estimate of the •, and F U is the phase unwrapping function. The first feature is the standard deviation of the IP, which is defined as IF may be easily derived using a simple backward finite difference (BFD) estimator of a form Similar to (6), the standard deviation of the IF has been considered as the next feature, such that with the difference that zero-centered version of the IF estimate is utilized to offset carrier frequency estimation errors, i.e., and scaling by the estimate of the frequency is applied to achieve invariance of signal bandwidth. BFD-based IF estimates suffer from their high variance. Hence, γ 3 feature is analogous to (8) albeit employs Kay's IF estimator with reduced variance [38] f (n) = 1 2π where After the initial experiments, υ was set to 63. Finite difference IF estimates of pulse radar signals are generally unreliable close to edges [24]. For this reason, features have been calculated after omitting 7% and 10% of the boundary samples for γ 2 and γ 3 , respectively. 2) IF Estimated by the QML Algorithm: In recent years, a number of more precise IF estimation methods have been developed in a field of time-frequency signal processing [39]. It can be observed that the NLFM can be accurately modeled as the polynomial-phase signal (PPS) [9], [30], [40]. This enables the application of the QML IF estimation method [41] tailored for the PPS class. In order to limit the computational burden, we skipped O'Shea refinement stage, so the QML method has been executed as follows.
It is assumed that s(t) is a PPS of a form where A denotes the constant amplitude, t ∈ [−T /2, T /2], and P is a degree of polynomial with the coefficients {a i |i ∈ [0, P ]}. In this research, we adopted the 12th degree PPS model [24], [31]. The IF of s(t) may be represented as Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply.
First, a series of discrete short-time Fourier transforms (STFTs) of x(n) using windows w h (m) of a length h from set H are obtained where DFT signifies discrete Fourier transform, m denotes discrete lag, k is the discrete frequency index, w h (m) = 1 for |m| < h/2, and w h (m) = 0 otherwise. During this study, the windows defined by {h ∈ H : h = 2 i ∧ h ≤ N }, i ∈ N were utilized. Next, IF estimates corresponding to consecutive STFTs are extracted from the ridge of the transform, i.e., Then, PPS coefficientsâ h = [â h 1 , . . . ,â h 12 ] T are computed for eachf h (n) by means of polynomial fitting. Since f (n) is normally odd function, it has been accounted that where y h , N h , and X h are defined in (17), at the bottom of the page.
T is the one which maximizes QML criterion of the formâ Finally, the IF of x(n) may be reconstructed fromâ f using (13). The γ 4 feature derived using the QML method is defined analogously to (8).

B. Moment of the Complex Envelope
Moments and cumulants of the complex envelope characterize distribution of a signal constellation. Hence, they are particularly useful for a recognition of phase-modulated waveforms. Yet, these metrics proved to be relevant also for modulations other than PSK [8].
The pth-order mixed moment of the complex random process with zero mean is given by [42] where E denotes the excepted value, * signifies the complex conjugate, and q is an integer smaller than p/2. Typically, the above quantity is estimated using the formula [2], [8] After the initial simulations, M 10 estimator that maximized interclass separation was chosen as γ 5 feature.

C. Features Derived From TFR
Among many existing time-frequency decompositions, the Wigner-Ville distribution (WVD) and the Choi-Williams distribution (CWD) are most prevalent in radar AMC. The continuous form of the WVD is given by (22) where τ denotes lag in continuous time domain. WVD is known to be optimal transformation for LFM analysis [43] in the sense that it follows its IF, i.e., where δ(f ) denotes the delta function. However, for higher order PPS with P > 2, WVD reveals inner cross-terms, which appears as artifacts in the midway of the chord linking two arbitrary points of the IF curve, as depicted in Fig. 5. This property provides the opportunity to distinguish between LFM and NLFM with TFR concentration measures.
To compute the feature, first discrete WVD needs to be derived [44] ρ wv (n, k) = then, L 4 /L 2 norm metric proposed by Jones and Parks [45] may be employed To retrieve the IF of higher order or multicomponent signals, the reduced interference distributions (RID) are preferable to WVD. The CWD, also known as exponential distribution, belongs to this class. As any of the Cohen class distributions, the CWD may be represented as a radar ambiguity function (RAF) multiplied by a kernel function Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply. where ν is the frequency lag, A s (ν, τ) denotes the RAF of s(t), and scaling factor σ controls the tradeoff between interference suppression and frequency resolution. During the research, we applied σ = 0.05, likewise to [2]. The discrete equivalent of (26) may be computed directly from the WVD of x(n) in a straightforward manner by executing 2-D DFT followed by 2-D inverse DFT where g cw (l, m) denotes the discrete CWD kernel and l is the discrete frequency lag. This operation has been illustrated in [24].
In order to extract the feature, the time-integrated CWD is required so that ρ f cw (k) kurtosis could be computed where K is a number of frequency-domain samples, andρ f cw and σ ρ * denote the mean and standard deviation of ρ f cw (k), respectively.

D. Spectral Features
Spectral features extracted from PSD are often employed in AMC. PSD estimatesP (k) are generally derived using the periodogramP Klauder et al. [46] proved that for a plain chirp signal with time-bandwidth product (TBP) exceeding 100, more than 99% of signal energy is confined between [−δF/2, δF/2] and its spectrum is asymptotically rectangular. This is not the case for the NLFM waveforms, which spectra are bell-shaped in order to attenuate the sidelobes. Spectral flatness known from acoustic signal processing [47] can be utilized to exploit this property .
The above metric takes values from 0 for a single harmonic to 1 for white noise. This feature is determined on a PSD portion corresponding to bandwidth B, as shown in Fig. 6.
Periodogram estimates have relatively high variance. To reduce the variance, alternative Welch's method was utilized for γ 9 . The PSD estimation was performed on eight segments with 50% overlap and windowed using Hamming function.

E. FrFT-Based Feature
In the past few years, the FrFT has drawn much attention in the field of radar signal processing [20], [48], [49]. The continuous form of the FrFT is defined by the formula [50] where u signifies transform domain, α ∈ R denotes transform angle, n i ∈ Z, and Intuitively appealing interpretation of the FrFT is the counterclockwise rotation of a signal in the time-frequency plane by an angle α. In case of the LFM waveform, there exists a single rotation angle α opt that results in energy integration in the transform domain [48]. This optimal angle is related to the chirp rate µ with the equation [20] α opt = 2 π arctan(µ).
As the result, the chirp is compressed into a single harmonic, as shown in Fig. 7. The same does not apply to the NLFM, where chirp rate µ(t) is a time-dependent function. Hence, effective compression occurs for a wider span of FrFT angles. We took advantage of this fact to distinguish between both waveforms. The feature extraction algorithm is provided in the following.
At first, compute a series of discrete FrFTs of the x(n) denoted as X ′ α (r), where r stands for discrete FrFT domain. Consecutive transforms are computed using rotation angles from the set where δ c α is the coarse search step. This α opt blind search strategy allows to circumvent the requirement for chirp rate estimation. Then, the maximum of magnitude of each FrFT is determined as what facilitates rough estimation of the optimal rotation angle using the obtained coarse FrFT profilê Next, by repeating the same procedure with the set A ′′ s.t. {α ∈ A ′′ : α = iδ f α ∧[α opt −(∆α/2) ≤ α ≤α opt +(∆α/2)]}, where δ f α means fine search step, i ∈ N, and ∆α is search span, new FrFT set X ′′ α (r) is obtained. Consequently, fine FrFT profile X ′′ max (α) is derived analogically to (35), so that α opt estimate might be improved.
Lastly, the feature is a skewness of the FrFT profile X ′′ max (α), which is calculated as and is interpreted in Fig. 8. It should be noted that X ′′ max (α) should be normalized by its maximum beforehand in order to make the feature SNR-invariant.
The definition (37) applies to upchirp waveforms, whereas the reverse is taken in case of downchirp. Following this rule, the feature value is always confined within an interval [0, 1]. During the simulations, the parameters δ c α = 1 • , δ f α = 0.1 • , and ∆α = 5 • were utilized, as a compromise between processing complexity and feature relevance.

V. FEATURE SELECTION
Once the features have been defined, the feature selection was performed. First, it enables the use of a simpler classifier, hence reduces the overall computational demand of the AMC method. Second, the curse of dimensionality which is inherent for a large set of features can be avoided due to this process.
Let us denote the original set of n Γ = 10 features as Ξ = {γ 1 , γ 2 , . . . , γ nΞ }. The problem of feature selection comes down to picking such a subset Γ * of a minimum possible cardinality n Γ < n Ξ , which maximizes some objective function J(·), i.e., The overall classification accuracy related with the subset Γ * is often comparable to a hypothetical accuracy with all available features employed. It is due to the fact that irrelevant and redundant variables are eliminated during the selection stage.
The feature selection was divided into two steps. During the first two, the features were evaluated individually against noise and multipath influence. This explorative analysis makes the results applicable for other scholars. The terminal step utilized minimal-redundancy-maximal-relevance (mRMR) criterion to choose the target subset Γ * for the designed algorithm.
The simulations presented in Sections V-VII were carried using a set of complex, baseband test signals buried in AWGN, with sampling frequency f s set to 100 MHz. The necessary code was executed in MATLAB 2021b. As a general rule, we utilized the proprietary software, except the embedded MATLAB functions periodogram and pwelch and the FrFT implementation developed by Ozaktas et al. [51].
To measure the relevance of the proposed features, their respective Fisher scores [52] were estimated. Let n c = 2 denote the total number of classes and y be the class index, 1 for LFM and 2 for NLFM. During each simulation, n x samples of x(n) were generated. For a specific ith feature, this yielded the dataset γ i ∈ R 1×nx . Consequently, the Fisher score of this individual feature has been approximated using the formula where n x|y means a number of samples belonging to the yth class,γ i and σ 2 i = nc y=1 n x|y σ 2 i|y are the mean and standard deviation of the ith feature estimated from γ i , whereasγ i|y together with σ 2 i|y are the same metrics obtained for the yth class only.
During the feature selection, n x = 1000 signal samples were utilized to compute (39) value, divided equally between LFM and NLFM C, i.e., n x|1 = n x|2 = 500. This applies to a given value of test parameter (single point in the graph). Fixed pulse duration T = 40 µs and bandwidth B = f s /40 were applied.
First, the noise influence on respective features has been examined. To this end, the Fisher score as a function of SNR was determined, in the range of −10 ÷ 30 dB with 1.25 dB steps. The results are shown in Fig. 9. Features derived from the complex envelope are marked with black, TFR-related features are plotted green, spectral features are represented with blue, QML feature is plotted in yellow, and the single FrFT feature is drawn using red color. It should be noted that  common temporal features, especially those derived from the IF, are highly susceptible to noise, what corroborates earlier studies, e.g., [24], [53].
The same experiment was conducted with R being a test parameter, so as to assess whether the features will keep their properties in multipath channels. The results are presented in Fig. 10. Not surprisingly, the features become less informative, while a number of reflections increases. One can infer that features derived from linear decompositions, such as STFT (γ 4 ) and FrFT, are less prone to interference in comparison to those extracted using bilinear transforms due to cross-terms inherent to the latter. This has also been noticed in [40] and [48]. Interestingly, WVD feature turned out to be more informative than one derived from CWD. This may be due to the fact that RID class distributions reduce not only crossterms but also signal autoterm, so can degrade overall output SNR [54]. It is also worth emphasizing that a couple of the abovementioned metrics are perceived as suitable for radar AMC, including PSD, BFD-estimated IF, or various TFRs combined with CNN classifiers. Meanwhile, they are multipath sensitive what may affect the overall accuracy of a designed AMC method. Fig. 11 gives an insight into this issue, where three different characteristics of a real-life NLFM waveform are compared, whereas one recording is barely distorted and the other is heavily corrupted due to the multipath effect.
During the preliminary selection, we considered the features which the score falls below 1 as insufficiently relevant. In other words, if F (γ) < 1 in a given range of test parameter, the feature was eliminated. First, we assumed that SNR ≥ 0 dB requirement is entirely sufficient for a typical ELINT/ES scenario, so that γ 1 , γ 2 , and γ 3 could be excluded. In the case of multipath, we wanted the features to keep the minimal Fisher score for R = 2, so γ 5 together with γ 7 was further discarded. Furthermore, among two spectral metrics γ 8 and γ 9 , less informative γ 8 was eliminated. Finally, γ 4 , γ 6 , γ 9 , and γ 10 remained for the last stage of feature selection.
The already mentioned mRMR algorithm [55] was chosen to determine the optimal subset of features. The mRMR not only favors the subsets, which increase the interclass separation, but also depreciates the combinations of similar features.
(40)  The relevance of a subset was estimated from Γ as the mean value of the mutual information I(·) between the samples related to consecutive features γ i and the dataset class label vector Ψ = [y 1 y 2 · · · y nx ] ∈ {1, 2}, i.e., Consequently, the redundance was estimated as mean mutual information between all existing feature pairs Given both values are known for all considered Γ, the mRMR-optimal subset if the one which maximizes the criterion (43) In this stage, fixed SNR = 10 dB and R = 5 were applied. To calculate the mutual information, we utilized functions from the repositories [56], [57]. The results for all existing subsets n Γ ≥ 2 of the remaining four features are shown in Fig. 12. The mRMR score is marked with red, and relevance together with redundance are plotted using dark and light gray, respectively. It can be noted that time-frequency and spectral metrics (WVD, PSD, and FrFT) are strongly correlated. At the same, the only temporal, IF-based γ 4 feature turned out to be the most discriminative.
According to the results, the optimal Γ * = {γ 4 , γ 9 } subset should be selected. However, γ 9 prerequisites bandwidth estimation, so would certainly deteriorate its relevance in practice. For this reason, the suboptimal {γ 4 , γ 10 } combination of QML and FrFT features was chosen as the ultimate set.

VI. WAVEFORM CLASSIFIER
Features scrutinization discovered the highly discriminative signal representation. This, in turn, reduced the requirements for the classifier's performance. In order to find the proper decision method, we synthesized sample signal set and compared outcomes of five popular classifiers, namely, linear and quadratic discriminants, support vector machine, k-nearest neighbors, and the NB. After the 10-fold cross-validation supervised test, it turned out that all of the above methods exhibited near perfect accuracy, except the inferior linear discriminant. Finally, Gaussian NB classifier was employed in the designed AMC algorithm due to its simplicity.
Given new observation of x(n) is associated with the true class label y and the extracted features are confined in multivariate random variable X = (γ 4 , γ 10 ) T , the classifier's task is to assign a such class-label y * , which will maximize the functional of a form [58] y * = arg max y∈{1,2} P (y|X).
Explicit solution of this problem would require the estimation of the prior probability P (y|X) what is infeasible, yet its equivalent form may be derived using the Bayes theorem P (y | X) = P (X | y) P (y) P (X) = P (γ 4 , γ 10 | y) P (y) P (γ 4 , γ 10 ) where P (X | y) is posterior and P (y) represents likelihood. The denominator in the above equation stays constant regardless of considered y, hence can be discarded. Furthermore, the NB classifier naively assumes that all features make equal and independent contribution to the result. Under these two conditions, the above problem can be further simplified to where ∝ signifies proportionality.
It can be noticed that if any of the conditional probabilities is close to zero, the arithmetic underflow may occur. In order to prevent this, the above formula is usually computed in log space The remaining problem is to estimate the conditional probabilities. Gaussian NB, as the name suggests, assumes that features follow the normal distribution, i.e., It can be deduced that NB supervised learning comes down to determine the parameters, which governs (48) distribution. This can be achieved using the training dataset along with sample mean and sample standard deviation estimators.
For this purpose, we synthesized n x = 10 000 training samples with n x|1 = 5000 and n x|2 = 5000, whereas the latter class was divided equally between NLFM A-D. Please note that in such case, P (y = 1) = P (y = 2), so that prior probability in (47) might be omitted. The remaining parameters of training dataset were as follows: SNR ∼ U (0, 20) dB and TBP ∼ U (50, 500) with B governed by U (f s /80,f s /8) and T bounded between [5 µs, 100 µs] at once. These values are typical of real-life reconnaissance scenarios. Finally, R = 0 was set during the training to avoid possible data corruption arising from impreciseness of multipath model (let us recall that multipath presence has been already accounted during the feature selection stage). The SNR, TBP, pulse duration, and bandwidths were modeled as random variables so as to reflect the fact that signals emitted by the noncooperative radars normally have unknown and varying parameters.
The scatter diagram representing sample dataset is depicted in Fig. 13. NLFM and LFM samples are plotted using green and red color, respectively. Black dashed line indicates obtained decision boundary. The samples related to respective NLFM waveforms are marked using various shapes.

VII. RESULTS AND DISCUSSION
In order to evaluate the algorithm, we applied two-pronged approach. First, classic quantitative research using softwaregenerated signals was conducted. This was followed by the experiment using real-life signals. The latter part acts as proof of concept, demonstrating that a properly designed AMC method is able to cope with multipath-distorted input signals.
For the sake of comparison, we modeled a part of real-time AMC method by Iglesias et al., tailored to field-programmable gate arrays (FPGAs). The NLFM/LFM pair is discriminated therein using the mean squared error between unwrapped, BFD-estimated IP and its least squares quadratic fitting. We did not modeled subblock processing, treating all x(n) samples as a whole. The decision threshold th LFM = 0.4906 was determined using training dataset with the same  parameters as provided in Section VI, with the except that SNR followed U (10, 20) dB to preserve BFD estimates from possible degradation. An interested reader can find necessary implementation details in [16].
The results are presented in the following using accuracy metric and confusion matrices. The accuracy is defined herein as where cŷ y denotes a total number of samples belonging to the yth class, which were assigned to theŷth class.

A. Simulation Results
This subsection shows the accuracy curves in a function of SNR and R. Single datapoint was computed using n x = 2000 test samples, divided equally between LFM and NLFM A-D, i.e., n x|1 = n x|2 = 1000. Again, TBP followed U (50, 500) distribution with B ∼ U (f s /80,f s /8) and T bounded between [5 µs, 100 µs] at once. Fig. 14 presents the accuracy of two considered algorithms as the SNR function, calculated with 1.25 dB steps and R = 0. Confusion matrix related to this experiment is given in Table II.
Please note that if a decision error is biased toward one class given equally distributed test data, the accuracy will never drop below 50%, even though an algorithm completely loses its discriminative properties.
Perhaps not surprisingly, the proposed method exhibits better noise resilience with 18-dB improvement when compared to ISP-based Iglesias method. This is obviously due to applying more complex features, but more importantly, the maximum accuracy of our method remains close to 100% for high SNR values, performing 5% better with respect to the ISP algorithm. Given the above results, next experiment was conducted with SNR confined between U (15, 25) dB so as to assure a fair comparison. In order to assess the multipath robustness, this time the accuracy was calculated as a function of R ∈ [0, 35] using the multipath model introduced in Section III. The results are given in Fig. 15 and Table III.
It is worth reminding that the proposed multipath model is not based on empirically determined channel parameters. As a consequence, the obtained results should be rather interpreted qualitatively. Nonetheless, it is clear that our algorithm is less prone to multipath and shows a potential to maintain stability in unfavorable propagation environments. This was verified in the next stage of research.

B. Results on Real-Life Signals Set
In the final step of the research, we did a measurement campaign and acquired signals originating from five different types of noncooperative, land-based pulse radars. This set of sensors included long/medium range surveillance and ATC radars, which utilized L and S frequency bands. Three of the emitters employed unknown NLFM waveform and the remaining two utilized linear chirp. Signals were intercepted from the ground with minor, natural and artificial obstacles present along the propagation path. This obviously introduced multipath components to received composite signal.
The data have been recorded using USRP-2901 software defined radio with omnidirectional antenna, controlled with LabView 2016. Baseband sampling frequency was set to 7.7 MHz. Each time after the data collection has been finished, samples were stored on disk drive for offline processing using MATLAB 2021b.
First, the necessary preprocessing was performed, starting with increasing the sampling rate with interpolation by factor 13. This adjusted f s to ca. 100 MHz-the value applied during simulation research, so that the same software could be used. Next, pulse detection was required. This was done  For each radar, a total of 100 pulses originating from antenna main lobe were selected, yielding the whole dataset consisting of 500 pulses. The essential processing consisted of running AMC algorithm with accordance to Sections IV and VI. Fig. 16 depicts four instances of recorded pulses. This exemplifies how prominent the multipath distortions could be. The instantaneous amplitude together with real and imaginary part of x(n) is plotted in gray, blue, and red, respectively.
Maximum, minimum, and median value of SNR for whole recorded dataset equaled −0.7, 53.2, and 21.1 dB, respectively, whereby SNR was estimated using the formula SNR ∼ = 10 log 10 (P s /σ 2 ε ) witĥ and the average noise powerσ 2 ε approximated as the mean squared value of 3N consecutive noise samples adjacent to the pulse. The confusion matrix estimated for the dataset is shown in Table IV.
Finally, the overall classification accuracy totaled 65.0% for Iglesias et al. method and 93.6% in case of our algorithm.

VIII. CONCLUSION
This concise introduced the novel, pattern recognition AMC algorithm for frequency-modulated radar waveforms. Unlike the existing methods, the algorithm is able to recognize NLFM signal with generalized frequency modulation law, so classification between LFM and NLFM depends not on the f (t) form but solely on the presence of nonlinear IF components. Designed method was successfully verified using real-life signals. The fact that both the QML IF estimator [59] and FrFT transform [49] have been already realized on FPGA makes the potential to run the proposed method in real-time systems.
We believe that more attention should be given to multipath propagation in radar AMC. While we acknowledge that the proposed phenomenon model may not be sufficient, it could still prove useful in prospective research. The alternative models of a multipath fading, such as Rician channel utilized in [15], should also be considered. It is worth emphasizing that currently celebrated CNN classifiers might be prone to multipath interference, depending on the class of applied TFR.
The obvious limitation of the algorithm is limited number of recognized waveforms. For further studies, we believe that it is particularly important to investigate modern, yet currently not recognizable low probability of intercept radar modulations, such as polyphase-coded FM (PCFM).