EP4334742A1 - Methods for computer estimation of a single tone and their application to radar systems - Google Patents

Abstract

Method for estimating the angular coordinates of point targets whose range has been estimated through a MIMO FMCW radar equipped with a plurality of transmitting (TX) and receiving (RX) antennas, wherein each couple of said TX and RX antennas is replaced with the equivalent virtual antenna; said virtual antennas form a virtual array composed by one or multiple horizontal uniform linear arrays (HULAs); said MIMO FMCW radar being arranged to generate real or complex signals in response to a propagation scenario including a plurality of point targets; said method including acquisition of a spectrum of said signal and its first derivatives, acquisition of a list of a predetermined number of discretized ranges, named frequency bins, of said previously detected point targets, and sequential execution for each of said frequency bins; said execution including: estimation of the horizontal spatial frequency and complex amplitude of the most dominant point target, the estimation including: computation of the spectrum its first three derivatives through EFT calculation, checking whether combination of the spectra associated with said Reference HULA, generates an amplitude peak in the resulting spectrum and when the checking is positive execution of the following steps: refinement of said range and complex amplitude of said most dominant point target; cancellation of said most dominant point target, calculation of a residual energy in the considered frequency bin and its comparison with a predetermined threshold, and when a residual energy exceeds is below the predetermined threshold computation of the spatial coordinates, otherwise updating of said acquired list of said frequency bins and return analyze a further frequency bin; finally, generation of an overall image of the propagation scenario.

Description

Methods for Computer Estimation of a Single Tone and Their Application to Radar Systems

Field of the invention

The present invention relates to a deterministic method for the detection of a tone and the estimation of its parameters in a colocated multiple-input multiple-output radar system operating in a slowly changing propagation scenario.

A tone, as wellknown, is a wave oscillation carachterized by a certain frenquency, amplitude and phase. Thus, the estimation of the single tone, means astimating its frequency, aplitude and phase.

State of the art

The problem of estimating the amplitude, phase and frequency of a tone in additive white Gaussian noise (AWGN) has received significant attention for a number of years because of its relevance in various fields, including radar systems and wireless communications. It is well known that the maximum- likelihood (ML) approach to this problem leads to a complicated nonlinear optimization problem. In practice, the most accurate ML-based single tone estimators available in the technical literature achieve approximate maximisation of this metric through a two-step procedure; the first step consists in a coarse search of tone frequency, whereas the second one in a fine estimation generating an estimate of the so called frequency residual (i.e., of the difference between the real frequency and its coarse estimate). Coarse estimation is always based on the maximization of the periodogram of the observed signal, whereas fine estimation can be accomplished in an open loop fashion or through an iterative procedure. On the one hand, all the open loop estimators exploit spectral interpolation to infer the frequency residual from the analysis of the fast Fourier transform (FFT) coefficients at the maxima of the associated periodogram and at frequencies adjacent to it as in "B. Quinn, "Estimation of frequency, amplitude, and phase from the DFT of a time series," IEEE Trans. Signal Process., vol. 45, no. 3, pp. 814-817, Mar. 1997." and "Q. Candan, "Analysis and Further Improvement of Fine Resolution Frequency Estimation Method From Three DFT Samples," IEEE Signal Process. Lett., vol. 20, no. 9, pp. 913-916, Sep. 2013.". Unlike iterative estimators, the accuracy they achieve is frequency dependent and gets smaller when the signal frequency approaches the center of one of the FFT bins. On the other hand, the iterative estimation techniques available in the technical literature are based on various methods, namely on: a) numerical methods for locating the global maximum of a function (e.g., the secan method "D. Rife and R. Boorstyn, "Single tone parameter estimation from discrete-time observations," IEEE Trans. Inf. Theory, vol. 20, no. 5, pp. 591- 598, Sep. 1974." or the Newton's method "J. Selva, "ML Estimation and Detection of Multiple Frequencies Through Periodogram Estimate Refinement," IEEE Signal Process. Lett., vol. 24, no. 3, pp. 249-253, Mar. 2017."); b) an iterative method for binary search, known as the dichotomous search of the periodogram peak as in "E. Aboutanios, "A Modified Dichotomous Search Frequency Estimator," IEEE Signal Process. Lett., vol. 11, no. 2, pp. 186- 188, Feb. 2004."; c) interpolation methods amenable to iterative implementation like "E. Aboutanios and B. Mulgrew, "Iterative Frequency Estimation by Interpolation on Fourier Coefficients," IEEE Trans. Signal Process., vol. 53, no. 4, pp. 1237-1242, Apr. 2005. " ; d) the combination of the above mentioned dichotomous search with various interpolation techniques as "Y. V. Zakharov, V. M. Baronkin, and T. C. Tozer, "DFT-based frequency estimators with narrow acquisition range," IEE Proc. - Commun., vol.148 , no. 1, pp. 1-7, Feb. 2001. e) the computation of the first derivative of the spectrum as in "C.-F. Huang, H.-P. Lu, and W.- H. Chieng, "Estimation of Single-Tone Signal Frequency with Special Reference to a Frequency-Modulated Continuous Wave System," Meas. Sci. Technol., vol. 23, no. 3, Mar. 2012.". Signal and System Models

In this session, we focus on the problem of estimating all the parameters of the complex sequence and its real counterpart with n = 0, 1, ..., N — 1 ; here, Ai and F₍ 6 0,1) denote the complex amplitude and the normalised frequency, respectively, of the l- th complex tone appearing in the right-hand side (RHS) of eq. (1), a₍ = |i4₍| , xpi = A{Ai) , w_{c n} is the n-th sample of an additive white Gaussian noise (AWGN) sequence (whose elements have zero mean and variance 2s² ) , w_{r n} = ¾{w_cn), N is the overall number of samples, and ¾{x) and arg(x) denote the real part and the phase, respectively, of the complex quantity x. It is useful to point out that the signal model (2) can be rewritten as x_r,n = åi=o [Qexp(/27mF_;) + Q^*exp(-y^'27mF₍)] + w_{r n}, ( 3 ) where represents the complex amplitude of the real tone appearing in the RHS of eq. (2).

The signal models (1) and (2) appear in a number of problems concerning radar systems, wireless communications and biomedical applications. In the remaining part of this paragraph, we analyse their meaning in the context of colocated and bistatic multiple- input multiple output (MIMO) radar systems operating in time division multiplexing (TDM) (see "S. M. Patole, M. Torlak, D. Wang, and M. Ali, "Automotive radars: A review of signal processing techniques," IEEE Signal Process. Mag., vol. 34, no. 2, pp. 22-35, Mar. 2017.") and at millimetre waves. These systems have the following important characteristics: a) They are equipped with a transmit (TX) antenna array and a receive (RX) antenna array, that consist of N_T elements and N_R elements, respectively; these allow to radiate orthogonal waveforms from different antennas and to receive distinct replicas of the electromagnetic echoes generated by multiple targets. b) Their TX and RX arrays are made of distinct antenna elements, placed at different positions. Moreover, TX antennas are close to the RX ones and, in particular, are usually placed on the same shield. c) Orthogonality of the transmitted waveforms is achieved by radiating them through distinct TX antennas over disjoint time intervals.

The architecture of the radar system considered in this work is illustrated in Fig. 1. In the remaining part of this section, two different models are described for the RF signal generated by the voltage controlled oscillator (VCO) of its transmitter and radiated by its TX array after power amplification. In the first case, corresponding to a frequency modulated continuous wave (FMCW) radar system, the VCO is fed by a periodic ramp generator; this produces a chirp FM signal, whose instantaneous frequency evolves periodically, as illustrated in Fig. 2. In this figure, the parameters T, T_R and T₀ represent the chirp interval, the reset time and the pulse period (or pulse repetition interval), respectively, whereas the parameters f₀ and B are the start frequency and the bandwidth, respectively, of the transmitted signal. If we assume that the radar system exploits all the available TX diversity (i.e., all its TX antennas) and that a time slot of T₀ s is assigned to each TX antenna, transmission from the whole TX array is accomplished over an interval lasting T_F = N_TT₀ s; this interval represents the duration of a single transmission frame.

Let us focus now on a single chirp interval and, in particular, on the time interval (0,G), and assume that, in that interval, the p-th TX antenna is employed by the considered radar system (with p 6 {0 , 1, ..., N_T — 1}); the signal radiated by that antenna can be expressed as

^S _RF(0 ⁼ A_RF ¾{s(t)}, (5) where A_RF is its amplitude, s(t)= exp(/0(t)), (6) and m = B/T is the chirp rate (i.e., the steepness of the generated frequency chirp). Let ^r _RP(.t^') denote the signal available at the output of the q-th RX antenna in the same time interval, with q = 0, 1, N_R — 1 (see Fig. 1). This signal feeds a low noise amplifier (LNA), whose output undergoes downconversion, filtering and analog-to-digital conversion at a frequency f_s = 1/T_s; here, T_s denotes the sampling period of the employed analog-to-digital converter (ADC). If the radiated signal s_RF(t) (5) is reflected by L static point targets, the useful component of (t) consists of the superposition of L echoes, each originating from a distinct target. In this case, if the propagation environment is static or undergoes slow variations, a simple mathematical model can be developed to represent the sequence of samples generated by the ADC in a single chirp interval. In deriving that model, the couple of the involved physical TX and RX antennas (namely, the p-th TX antenna and the q-th RX antenna) of the considered bistatic radar is usually replaced by a single virtual antenna of an equivalent monostatic radar. In particular, if it is assumed that all the TX and RX antennas belong to the same planar shield and a reference system lying on the physical antenna array is defined, the abscissa x_v and the ordinate y_v of the 17-th virtual antenna element associated with the p-th TX antenna and the q-th RX antenna can be computed as (e.g., see "J. Li and P. Stoica, Eds., MIMO Radar Signal Processing. Hoboken, NJ: J. Wiley & Sons, 2009. ") and

_ yp+yg yp ₂ (9) respectively, with v = 0, 1, N_VR — 1; here, are the coordinates of the TX (RX antenna) and N_VR=N_T-N_R represents the overall size of the resulting virtual array. Based on these assumptions, the n-th received signal sample acquired through the 17-th virtual antenna element (with v = 0, 1, ..., N_VR — 1) can be expressed as (e.g., see "J. Gamba, Radar Signal Processing for Autonomous Driving. Springer Singapore, 2020.") with n = 0, 1, ..., N 1; here, N is the overall number of samples acquired over a chirp period, a₍ is the amplitude of the Z-th component of the useful signal, is the normalized version of the frequency characterizing the Z-th target detected on the v-th virtual RX antenna, is the delay of the echo generated by the Z-th target and observed on the 17-th virtual channel, R_\, 0; and f_i denote the range of the Z-th target, its azimuth and its elevation, respectively, f\ⁿ⁾ = 2nf_Qr\^v) , ( 14 )

(nL and w_rn, the n-th sample of the AWGN sequence affecting the received signal, is a Gaussian random variable having zero mean and variance s² (assumed to be independent of v). It is important to point out that: a) the real signal model (10) is identical to that expressed by eq. (2) and can be adopted in all the FMCW radar systems that extract only the in-phase component of the signal captured by each RX antenna; b) some commercial MIMO radar devices provide both the in-phase and quadrature components of the received RF signals (e.g., see "M. A. Richards, Fundamentals of Radar Signal Processing. New York: McGraw-Hill, 2005."). In the last case, the complex model equivalent to that expressed by eq. (1), must be adopted in place of its real counterpart (10) for any n; here, for any v and Z.

The second case we consider in the generation of the radiated waveform corresponds to a stepped frequency continuous wave (SFCW) radar system represented in Fig. 3. Its name is motivated by the fact that the VCO of its transmitter is fed by a staircase generator. For this reason, the instantaneous frequency of the resulting RF signal takes on N distinct and uniformly spaced values in an interval lasting T s for each TX antenna (see Fig. 4); the n-th value of the instantaneous frequency is f_n=fo+nAf, (17) with n = 0, 1, , N — 1; here, f₀ is the minimum radiated frequency, D/ is the frequency step size and N is the overall number of transmitted frequencies. It can be shown that, under the same assumptions made in the derivation of eq. (15), the measurement acquired through the v-th virtual element at the n- th frequency can be expressed as with v = 0, 1, ..., N_{V R} — 1; here, the amplitude ip_t and the phase

_(V)

A_t are still expressed by eqs. (14) and (16), respectively, is the normalised delay characterizing the Z-th target and observed on the v-th virtual antenna, and the parameters (a_;, have exactly the same meaning as the one illustrated for the received signal model (15).

Then, in both the considered FMCW and SFCW radar systems, the signal observed on each virtual channel can be represented a superposition of L real or complex oscillations; moreover, the value of the parameter L has to be considered unknown. In the following derivations, the real samples n = 0, 1, ..., N — 1} or their complex counterpart n = 0, 1, ..., N — 1} acquired on the l?-th virtual channel are collected in the iV-dimensional vector with z = r or c. This vector is processed by the next stages of the radar receiver for target detection and estimation. As it can be easily inferred from eqs. (10)—(12) ((15)— (16) and (18)— (19)), in the considered radar system, the problem of target detection and range estimation on the v-th virtual channel is equivalent to the classic problem of detecting multiple overlapped sinusoids (multiple overlapped complex exponentials) in the presence of AWGN and estimating their frequencies. In fact, if, in a FMCW radar system, the Z-th tone is found at the frequency f_t , the presence of a target at the range (see eqs. (12) and (13)) p — ^lc f(^v)

^Kv ( ·¹ ~ 2 21) ~_m h is detected. Similarly, in a SFCW radar system, the normalised

—f_t?) delay F_; estimated on the the v -th virtual channel is associated with a target whose range is (see eqs. (13) and (19))

^Kp _ hd p( v,l ^v) 1 (22)

2 —Af G, .

Information about the angular coordinates (namely, the azimuth and the elevation) of this target, instead, can be acquired

_(V\ through the estimation of the set of N_VR phases {ip_t ; v = 0, 1, ..., N_VR — 1} observed over the available virtual antennas. In fact, since (see eqs. (13) and (14)) where is the wavelength associated with the frequency /₀, the sequence exhibits a periodic behavior characterized by the normalised horizontal spatial frequency if the considered virtual elements form an horizontal uniform linear array (ULA), whose adjacent elements are spaced d_H m apart. Dually, if a virtual vertical ULA is assumed, the periodic variations observed in the same sequence of phases are characterized by the normalised vertical spatial frequency

F_v,i A2^sin(0_;), (26) where d_v denotes the distance between adjacent elements of the virtual array itself. Consequently, angle finding can be easily accomplished by digital beamforming, i.e. by performing FFT processing on the estimated phases taken across multiple elements of the virtual array in a single frame interval.

From the considerations illustrated above, the following conclusions can be easily inferred: a) on the one hand, achieving precise estimation of the range of a given target requires the availability of an accurate estimate of the normalised frequency (or delay) of the real or complex tone associated with the target itself over at least one RX antenna; b) on the other hand, the quality of the estimate of the direction of arrival (DOA) characterizing the radar echo generated by a given target depends on the accuracy of the phase estimated over multiple virtual channels. The last consideration motivates the importance of accurately estimating this parameter over multiple antennas. Finally, in analysing the suitability of multiple tone estimators to colocated MIMO radar systems operating at millimetre waves, the following additional technical issues need to be taken carefully into account:

1) These radar systems often operate at short ranges and in the presence of extended targets. Each of resulting radar images is a cloud of point targets whose mutual spacing can be very small. For this reason, the accuracy of these images depends on the frequency resolution (delay resolution) achieved on each virtual antenna in a FMCW (SFCW) radar system. In fact, this makes the receiver able to separate point targets characterized by similar ranges.

2) Information about the RCS of each point targets can be exploited in the classification of extended targets; for this reason, the availability of an accurate estimate of the amplitude of each radar echo can be very useful in a number of applications (e.g., in SAR imaging).

3) Distinct radar echoes can be characterized by substantially different signal-to-noise ratios (SNRs).

4) The number N of samples acquired over each virtual channel usually ranges from few hundreds to few thousands.

The last two issues explain why significant attention must be paid to the accuracy achieved by the adopted estimation algorithms at low SNRs and/or for relatively small values of N, since this can appreciably influence the quality of the generated radar image.

Summary of the invention

Scope of the present invention is to provide for methods for estimation of a single tone, which has a low computation complexity and that can be easily implemented in low computation power such has radars and in the base stations of cellular communication networks.

The basic idea of the present invention is based on the approximate maximum likelihood estimation of the single tone. The present invention proposes a new method that can be declined into different ways according to the reference domain, such as real or complex frequency or time domain. This is the reason for which the same inventive concept is described with the aid of the claims 1 and 10. More in detail, when the single tone represents a frequency, the present invention proposes a method for computer estimation of the frequency, amplitude and phase of the single tone defining a signal acquired by means of a receiver, wherein the tone is real or complex; the method including

- estimation of a spectrum (X₀) of said signal and its first N_FFT — 1 derivatives ({X_fc;k = 1,...,N_FFT — 1}) with N_FFT = 3 when said signal is real and N_FFT = 4 when said signal is complex, through FFT calculation with an oversampling factor (M ),

- iterative calculation of an estimate of the parameters of the most dominant tone, including phase (y), amplitude (a) and frequency (/) on the basis of said spectrum, and its derivatives, through N_FFT — 1 FFTs.

When, instead, the single tone represents a delay, the present invention proposes a method for computer estimation of the delay, amplitude and phase of a single complex tone defining a signal acquired by means of a receiver; the method including

- estimation of the impulse response (X₀) characterizing the

— (v) communication channel and its first N_FFT — 1 derivatives ({X_fc ;k = 1,...,N_FFT — 1}) with N_FFT = 4, through IFFT calculation, with an oversampling factor (M),

- iterative calculation of an estimate of the parameters of the most dominant tone, including phase (ip), amplitude (a) and delay (t) on the basis of said impulse response, and its derivatives, through N_FFT — 1 IFFTs.

It is immediately apparent that the two implementations depicted above are dual concepts.

It is also clear that the receiver can be part of a radar or of any device capable to acquire a single tone signal. Thus, the present invention can be implemented also in 5G base stations. The present invetion aims at providing a total of three approaches about the estimation of a tone and their exploitation in colocated frequency modulated continuous wave (FMCW) and stepped frequency continuous wave (SFCW) multiple-input multiple-output (MIMO) radar systems operating at millimeter waves. In the present invention, a novel ML-based iterative estimator of a single real and complex tone is presented. The derivation of this estimator is based on: a) expressing the dependence of the ML metric on the frequency residual in an approximate polynomial form through standard approximations of trigonometric functions; b) exploiting the alternating minimization technique for the maximization of this metric. Moreover, its most relevant feature is represented by the fact that it requires the evaluation of spectral coefficients that are not exploited by all the other related estimation methods available in the technical literature.

It should be clear that the FFT represents the practical implementation of the DFT through computer means. Therefore, DFT and FFT can be used and referenced in an interchangeable way. The same concept applies also to IFFT and IDFT.

The present invention will be made fully clear with the aid of the attached dependent claim describing preferred embodiment of the invention.

Brief description of the figures

Figure 1 shows an example of architecture of a MIMO FMCW radar system; Figure 2 discloses the time diagram of the instantaneous frequency of the RF signal radiated in a FMCW radar system;

Figure 3 shows an example of architecture of a MIMO SFCW radar system; Figure 4 discloses the time evolution of the instantaneous frequency of the RF signal radiated in a SFCW radar system;

Figure 5 shows a block diagram of the SFE#1 algorithm;

Figure 6 shows a block diagram of the SFE#2 algorithm;

Figure 7 shows a block diagram of the CSFE#1 algorithm; Figure 8 shows a block diagram of the CSFE#2 algorithm;

Figure 9 shows a block diagram of the CSDE#1 algorithm;

Figure 10 shows a block diagram of the CSDE#2 algorithm;

Detailed description of exemplary embodiments of the invention

Approximate Maximum Likelihood Estimation of a Single Tone In this section, we derive a new method for estimating the parameters of a single real or complex tone.

Estimation of a single frequency

Let us focus first on the problem of estimating the parameters (namely, the frequency and complex amplitude) of a single tone contained in the real sequence {x_r,n; n = 0, 1, ..., N — 1}, whose n- th sample is expressed by eq. (10) with L = 1, i.e. as or, equivalently, as (see eq. (3)) where (see eq. (4)) and F are the complex amplitude and the normalised frequency, respectively, of the tone itself. It is well known that the ML estimates F_ML and C_ML of the parameters F and C, respectively, represent the solution of the least square problem

(F_ML< CML) = argmine(R,C), ( 30 )

F,C where F and C represent trial values of F and C_r respectively, is the mean square error (MSE) evaluated over the whole observation interval, is the square error between the noisy sample x_rn (28) and its useful component s_n(F,C)= Cexp (33) evaluated under the assumption that F = F and C = C.

The cost function e(R, C) (31) can be easily reformulated in a different way as follows. To begin, substituting the RHS of eq. (33) in that of eq. (32) produces, after some manipulation, f_h = 2phR, (35) and ¾{x} (3{x}) denotes the real (imaginary) part of the complex variable x . Then, substituting the RHS of eq. (34) in that of eq. (31) yields where

X_R{P)± M{X (F)}, X,(P)A 3{X(P)}, is, up to the scale factor 1/iV, the Fourier transform of the sequence {x_r,n},

Based on eq. (36), it is not difficult to show that the optimization problem expressed by eq. (30) does not admit a closed form solution because of the nonlinear dependence of the function e(R,C) on its variable P . However, an approximate solution to this problem can be derived by a) Exploiting an iterative method, known as alternating minimization (AM) (e.g., see "I. Ziskind and M. Wax, "Maximum likelihood localization of multiple sources by alternating projection," IEEE Trans. Acoust. Speech Signal Process., vol. 36, no. 10, pp. 1553-1560, Oct. 1988."). This allows us to transform the two-dimensional (2D) optimization problem expressed by eq. (30) into a couple of interconnected one^¬ dimensional (ID) problems, one involving the variable P only (conditioned on the knowledge of C), the other one involving the variable C only (conditioned on the knowledge of P). b) Expressing the dependence of the function a(F,C) on the variable F through the couple (F_c,5) such that

P = F_c + dF_dft, (40) where F_c is a given coarse estimate of F, d is a real variable called residual, is the normalized fundamental frequency associated with the N₀- th order discrete Fourier transform (DFT) of the zero padded version of the vector collecting all the elements of the sequence {x_r,n}, M is a positive integer (dubbed oversampling factor), 0_D is a D —dimensional (column) null vector and iV₀ = M ·N . c) Expressing the dependence of the function a(F,C) (36) on the variable d through its powers {dⁱ; 0 < l< 3}; this result is achieved by approximating various trigonometric functions appearing in the expression of s(F,C) with their Taylor expansions truncated to a proper order.

Let us show now how these principles can be put into practice. First of all, the exploitation of the above mentioned AM approach requires solving the following two sub-problems: PI) minimizing the cost function a(F,C) (36) with respect to C, given F = F; P2) minimizing the same function with respect to F, given C = C . Sub-problem PI can be easily solved thanks to the polynomial dependence of the cost function s(F,C) on the variables C_R and C_j . In fact, the minimum of the function s(F,C) with respect to the last variables is the solution of the equations and C,- X,(F)- C,g_R(P)+ C_R g,(F)= 0, (46) that result from computing the partial derivative of the RHS of eq. (36) with respect to C_R and C_Ir respectively. Solving the linear system of equations (45)-(46) in the unknowns C_R and C_j produces the optimal values and f-- -XR (F)gi(F)+Xl(F) [1+9R(F)]

(48) i-ls(F)l² of C_R and C_j , respectively. Putting together the last two formulas yields

Therefore, given F = F, the optimal value C of the variable C can be computed exactly through the last equation; this requires the evaluation of X(F) and g(F) (see eqs. (38) and (39), respectively). Note that, on the one hand, g(F) can be easily evaluated through its exact expression

1exp(-;^'4reNF)-l 5(0 7 (50)

N7 exp(-;^'4reF)-l’ which is easily derived from its definition (39). On the other hand, 7(F) can be computed exactly through its expression (38) or, in an approximate fashion, through a computationally efficient procedure based on the fact that the vector

X_s = M X₀ (51) collects iV₀ uniformly spaced samples of the function X(F) (38). In fact, the k-th element of the vector X₀ (42) can be expressed as where

F_fc= k F_DFT (54) with k = 0 , 1 , ..., iV₀ — 1 . For this reason, an approximate evaluation of the quantity X(F) at a normalised frequency F different from any multiple of F_DFT (41) can be accomplished by interpolating the elements of the vector X_s (51); the last vector, in turn, can be easily computed after evaluating the N₀- th fast Fourier transform (FFT) of C_O,ZR (43), i.e. the vector X₀ (42).

Let us take into consideration now sub-problem P2 . Unluckily, this sub-problem, unlike the previous one, does not admit a closed form solution. For this reason, an approximate solution is developed below. Such a solution is based on representing the normalized frequency F in the same form as F (see eq. (40)), i.e. as

F = F_c + d F_DFT (55) and on a novel method for estimating the real residual d, i.e. for accomplishing the fine estimation of F . This method is derived as follows. Representing the trial normalized frequency F according to eq. (40) allows us to express the variable f_h (35) as

Ft_i ^@ _n FtiA, (56) where

D = 2pd F_Q PT ( 57 ) is a new variable and

§_n = 2 ph F_c. (58)

Then, substituting the RHS of eq. (56) in that of eq. (34) (with C = C) yields that can be rewritten as

^£ _n(F, C) = X² _in + 2 [c² + Cf] +2(C| — C²)[cos(20_n)cos(2nA) — sin(20_n)sin(2nA)]

— 4C_RC_/[sin(20_n)cos(2n7\) + cos(20_n)sin(2nA)]. ( 60)

If the normalized frequency F_DFT (41) is small enough (i.e., if the FFT order N₀ is large enough) , the trigonometric functions cos (knK) and sin (knK) appearing in the RHS of eq. (60) (with k = 1 and 2) can be approximated as sin(nA) — nA, (61) and

Substituting the RHSs of eqs . (61) -(64) in that of eq. (60) and, then, substituting the resulting approximate expression in the RHS of eq. (31) yields

^£SFE(A C) = here, p = F_C/F_dft, 66 g_p\n\= n^p for p = 1, 2 and 3, g₀\n\= 1, for any p and k = 1,2, and

^Xk,n — 71 ·Cg-gi (69) with n = 0, 1, ..., N — 1. It is important to point out that: a) if p is an integer, the quantity X_kp (68) (with k = 1 and 2 ) represents the p-th element of the vector generated by the N₀-th order DFT of the zero padded version of the vector b) if p is not an integer, the quantity X_kp can be evaluated exactly on the basis of eq. (68) or, in an approximate fashion, by interpolating a subset of adjacent elements of the vector X_fc (70); c) the evaluation of the vectors {X_fc; k = 1,2} requires two additional FFTs.

Since the function e_5RE(D,C) (65) is a polynomial of degree three in the variable D, an estimate D of D and, consequently, an estimate (see eq. (57)) of d, can be obtained by computing the derivative of this function with respect to D, setting it to zero and solving the resulting quadratic equation a(r)D² + b(p)D+ c(p)= 0, (74) in the variable D; here, and

Note that only one of the two solutions of eq. (74), namely b(p)+J(b(p))^Z®a(p)c(p)

D (78) 2a(p) has to be employed. A simpler estimate of D is obtained neglecting the contribution of the first term in the left-hand side of eq. (74), i.e. setting a(p)= 0. This leads to a first- degree equation, whose solution is

D = —c{p)/b{p). (79)

Given an estimate D of D (and, consequently, and estimate d of d; see eq. (73)), the fine estimate of F can be evaluated on the basis of eq. (55).

The mathematical results derived above allow us to derive novel estimation algorithms. Our first estimation algorithm, called single frequency estimator #1, (SFE #1), for iteratively estimating the normalized frequency F and the complex amplitude C is initialised by

1) Evaluating: a) the vector X₀ (42); b) the initial coarse estimate F^ of F as

F_c ⁽⁰⁾= aF_DFT, (81) where the integer ά is computed by means of the well known periodogram method (e.g., see "E. Aboutanios and B. Mulgrew, "Iterative Frequency Estimation by Interpolation on Fourier Coefficients," IEEE Trans. Signal Process., vol. 53, no. 4, pp. 1237-1242, Apr. 2005."), i.e. as ά = are c) the initial estimate C® of C on the basis of eq. (49) with F = ; d) the spectral coefficients and X₂p ^on the basis of eq. (68); e) the coefficients {K_p(2a); p = 1, 2} ({K_p(2a ); p = 1, 2, 3}) according to eq. (67), the coefficients (a(a), ό(a), c(a)} ( {&(<$) , c(a)} according to eqs. (75)-(77) and the first estimate D^⁰-¹ of D on the basis of eq. (78) (eq. (79)); f) the first fine estimate F® of F on the basis of eqs. (55) and (73), i.e. as

2) Setting its iteration index i to 1.

Then, an iterative procedure is started. The i-th iteration is fed by the estimates of F and C, respectively, and produces the new estimates F® and of the same quantities (with i = 1 , 2 , . .. , iV_SFE , where iV_SFE is the overall number of iterations); the procedure employed for the evaluation of F® and C® consists of the two steps described below (the p-th step is denoted SFE#1-Sp). SFE#1-Sl - The new estimate D® of D is computed through eq. (78) (eq. (79)); in the evaluation of the coefficients (a(p), b(p), c(p)} ({b(p) , c(p)}) appearing in the RHS of these equations, is assumed. Then, is evaluated.

SFE#1-S2 - The new estimate C® of C is evaluated through eq. (49); F = F® is assumed in this case. Moreover, the index i is incremented by one before starting the next iteration.

At the end of the last (i.e., of the iV_SFE-th) iteration, the fine estimates F = F(^WSFE) and C = C^^Nsfe^ of F and C, respectively, become available.

It is important to point out that: a) The estimate 5® of the residual d computed by the SFE#1 in its i-th iteration is expected to become smaller as i increases, since F® should progressively approach F if our algorithm converges. b) The coefficients {K_p(2a); p = 1, 2, 3} can be computed exactly on the basis of eq. (67). However, since the definition (67) can be put in the equivalent form where the identities and

+(3N³ - 3 N² -3 N - 1 )q^N+1 (89)

+( —37V³ + 6 N² - 4 )q^N+2 + (N - 1 )³q^N+3 (90) holding for any q 6 (C, can be exploited for the efficient computation of K_v(2d) for any p and a. c) The quantities k = 1,2} required in the first step of the i-th iteration can be computed exactly on the basis of eq. (68). However, they can be also evaluated, in an approximate fashion, by interpolating / adjacent elements of the iV₀ - dimensional vectors X_fc (70), where / denotes the selected interpolation order. d) The estimate D® evaluated according to eq. (79) is expected to be less accurate than that computed on the basis of eq. (78). However, our numerical results have evidenced that both solutions achieve similar accuracy. e) The SFE#1 can be employed even if the single tone appearing in the RHS of eq. (27) is replaced by the superposition of L distinct tones (see eq. (3)). In this case, the strongest (i.e., the dominant) tone is detected through the periodogram method (see eq. (82)) and the parameters of this tone are estimated in the presence of both Gaussian noise and the interference due to the remaining tones. Therefore, the estimation accuracy of the SFE#1 is affected by both the amplitudes and the frequencies of the other (L— 1) tones. f) A stopping criterion, based on the trend of the sequence {D® ; /= 1, 2, ...}, can be easily formulated for the SFE#1. For instance, the execution of its two steps can be stopped if, at the end of the i-th iteration, the condition

|D®|<e_D (91) is satisfied; here, e_D represents a proper threshold.

Our second estimation algorithm (called single frequency estimator #2, SFE#2) relies on the same mathematical results exploited in the derivation of the SFE#1 and has the same structure. However, it is also based on the idea of taking the estimate generated in its (i— l)-th iteration as a coarse estimate of F in the following (i.e., in the i-th) iteration. This leads to representing the (unknown) normalised frequency F as

F = F^(i_1) + 5F_dft, (92) where d is a residual whose estimate 5®, generated in the i-th iteration of SFE#2, is expected to tend to zero as / increases. The initialization phase of SFE#2 is very similar to that of SFE#1, the only difference being represented by the fact that the quantity is also defined. The first step of SFE#2 is different from the corresponding step of SFE#1 and can be summarised as follows.

SFE#2-Sl) The new estimate D® of D is computed through eq. (78) (eq. (79)); in the evaluation of the coefficients (a(p), b(p), c(p)} ({b(p) , c(p)}) appearing in the RHS of these equations, C = CV-V and are assumed. Then,

F⁽ⁱ⁾ = F⁽ⁱ-ⁱ⁾ +— (95)

2TC is evaluated.

The second step of SFE#2, instead, is the same as that of SFE#1.

Preferably, the method for the estimation of the frequency, the amplitude and the phase of a real sequence (see Steps 2.1.0 - 2.1.2 of Fig. 5) can be described as follows:

- (Step 2.1.0) representation of a tone normalized frequency F = f/f_{s r} where f_s is the sampling frequency, as sum of a main portion (QCF_dft) and a residual portion (5F_dft), where F_DFT = l/iV₀, where N₀ is the FFT order, wherein said main portion is determined by searching the index (a) of the vector element X_oa corresponding to the maximum absolute value in a first half plus one elements of the same vector X₀ and said residual frequency is found through an iterative procedure together with the corresponding tone complex amplitude C= aexp(jxp)/2, where a is said tone amplitude and y is said tone phase.

- (Step 2.1.1) said iterative procedure includes a preliminary initialization comprising: a) setting the iteration index i to 1, the initial estimate of tone complex amplitude C equal to said vector element (X_oa)_/ the initial estimate of A = 2pdF_DFT to zero and the initial estimate of F to said main portion (CCF_dft); b) computing the complex coefficients (K_p(2a);p = 1,2,3} as the 2a- th element of the iV₀-order FFT of the sequence n^p, with p = 1,2,3; then iteratively

- (Step 2.1.2) the i-th iteration is fed by the estimates of D, F and C computed in the previous iteration, to generate new estimates of D , F and C , for a predetermined number of iterations, wherein each iteration includes a) a first step where the residual D is calculated solving the first order equation Z¾D + c¾ = 0 OR the second order equation a_dD² + ί¾D + C_Q ⁼ 0 in the variable D with: where ¾{x} denotes the real part and 3{x} denotes the imaginary part of the complex quantity x; discarding one of the two solutions, which does not satisfy the inequalities: if a = 0 otherwise if said second order equation is solved; then the new estimate of F is obtained by adding said residual (A F_dft) to said main portion (af_DFT); b) a second step wherein said new estimate of F is employed to compute the new estimate of C as where is generated by interpolating the elements of said vector MX₀ where M is said oversampling factor employed in said FFT calculations.

The above mentioned Steps 2.1.1 - 2.1.2 of Fig. 5 can be replaced by an alternative method including (see Steps 2.2.1 - 2.2.2 of Fig. 6):

- (Step 2.2.1) said iterative procedure includes a preliminary initialization comprising: a) setting the iteration index i to 1, the initial estimate of C equal to said vector element (^o_,a)_/ the initial estimate of A = 27T5 F_DFT to zero and the initial estimate of F to said main portion (aF_DFT); b) computing the complex coefficients (K_p(2a); p = 1,2,3} as the 2a- th element of the iV₀-order FFT of the sequence n^p, with p = 1, 2, 3 ; c) computing the real quantity p = F/F_DFT; then iteratively - (Step 2.2.2) the i-th iteration is fed by the estimates of D, F and C and the real quantity p computed in the previous iteration, to generate new estimates of D, F, C and p , for a predetermined number of iterations, wherein each iteration includes a) a first step where the residual D is calculated solving the first order equation Z¾D + c¾ = 0 OR the second order equation a_dD² + Z¾D + C_Q ⁼ 0 in the variable D with: where ¾{x} denotes the real part, 3{L:} denotes the imaginary part of the complex quan e computed as: with p = 1, 2, 3 and i^x _k,p> ^k = 1-2} as: where x_rn is the n-th element of the said real signal; said coefficients ({X_kp;k = 1,2}) can be also computed by interpolating the elements of said spectrum (X_fc); discarding one of the two solutions, which does not satisfy the inequalities:

(0 < D < 7 TF_DFT if a = 0

—TTF_QPT < D < TTF_QPT otherwise if said second order equation is solved; then the new estimate of F is obtained by adding said residual (A F_dft ) to said estimate of F; b) a second step wherein the new estimate of p is computed by dividing said new estimate of F by F_dft ; c) a third step wherein said new estimate of F is employed to compute the new estimate of C as: where is generated by interpolating the elements of said vector MX₀ where M is said oversampling factor employed in said FFT calculations.

All the results developed above refer to the real sequence {x_r,n}, whose n-th element is expressed by eq. (27). However, similar estimation method - dubbed complex SFEs, CSFEs - can be developed for the complex counterpart, i.e. for a complex sequence {x_c,n} such that (see eq. (15) with L = 1 ) x_c,n = Aexp(j2nnF) + w_{c n}, (96) with n = 0, 1 , N — l . In this case, the ML estimation of the parameters F and A can be formulated in a similar way as eq. (30), the only differences being represented by the fact that: a) the parameter C is replaced by A ; b) the term s_n(F, C) appearing in the RHS of eq. (31) is replaced by where represents the useful component of x_{c n} (96) under the assumption that F = F and A = A . Substituting the RHS of the last equation in that of eq. (97) and, then, the resulting expression in the RHS of the MSE (see eq. (1)) yields, after some manipulation, where

The procedure adopted to minimise the last function with respect to F and A is similar to that illustrated above for the SFE. For this reason, in the following we limit to show essential mathematical results only. First all, it is easy to show that, given F = F, the function s{F,A) (100) is minimised with respect to A selecting

A=A = X(F), (102) where X(F) can be computed exactly through its expression (see eq. (38)) replacing the real sequence {x_r>n} with the complex sequence {x_c,n} ) or, i^{n an} approximate fashion, through an interpolation technique following the same mathematical approach as that illustrated for the SFE. Note also that the mathematical result expressed by eq. (102) is exact. On the contrary, given A = A, a closed form expression for the value of F minimising the function s(F,i4) cannot be derived because of the nonlinear dependence of this function on F. However, following the same mathematical approach as that illustrated for the SFE, the approximate expression can be obtained for the cost function s{F,A) (100); the quantities appearing in the last equation have the same meaning as the ones defined for eq. (65). It worth stressing that: a) The considerations expressed about the evaluation of the quantities {X_k,p} in our derivation of the SFEs apply to the CFSEs too. However, in this case, three additional spectral coefficients, namely {X_{k p}; k = 1,2,3}, need to be computed. b) The initial estimate p^ = a (see eq. (93)) of p is evaluated in a similar way as the SFEs (see eq. (82)), i.e., as a ^{= ar}g„ <ZE{0m,1,.a..,xN_Q —1J}*o,a|· (¹⁰⁴)

The minimization of the function ^CSFEC^W⁴) (103) with respect to D is achieved by taking its partial derivative with respect to this variable and setting it to zero; this results again in the quadratic equation (74), whose coefficients are b{p)= ¾{i4^*X_2,p} (106) and c(p)A -3{i¾ _,p}. (107)

Following the approach adopted in the development of the SFEs and exploiting the mathematical results expressed by eq. (102) and by eqs. (105)-(107) allow us to easily derive the CSFEs. Preferably, the method for the estimation of the frequency, the amplitude and the phase of a complex sequence (see Steps 2.3.0 - 2.3.2 of Fig. 7) can be described as follows:

- (Step 2.3.0) representation of a tone normalized frequency F = f/f_{s r} where f_s is the sampling frequency, as sum of a main portion (OCF_dft) and a residual portion (5F_dft), where F_DFT = l/iV₀, where iV₀ is the FFT order, wherein said main portion is determined by searching the index (a) of the vector element X_0Q corresponding to the maximum absolute value of the same vector X₀ and said residual frequency is found through an iterative procedure together with the corresponding tone complex amplitude A = aexp(jxp), where a is said target amplitude and y is said target phase.

- (Step 2.3.1) said iterative procedure includes a preliminary initialization comprising: setting the iteration index i to 1, the initial estimate of A equal to said vector element (Xo_,a)_t the initial estimate of A = 2i F_DFT to zero and the initial estimate of F to said main portion (CCF_dft); then iteratively

- (Step 2.3.2) the i-th iteration is fed by the estimates of D, F and A computed in the previous iteration, to generate new estimates of D , F and A , for a predetermined number of iterations, wherein each iteration includes a) a first step where the residual D is calculated solving the first order equation Z¾D + % = 0 OR the second order equation in the variable D with: a_a = bji4 X_3Qj, where ¾{x} denotes the real part and 3{x} denotes the imaginary part of the complex quantity x; discarding one of the two solutions, which does not satisfy the inequalities: if a = 0 otherwise if said second order equation is solved; then the new estimate of F is obtained by adding said residual (A F_dft) to said main portion (af_DFT); b) a second step wherein said new estimate of F is employed to compute the new estimate of A as A = Xi_nt{P) and Xi_ntif^) is generated by interpolating the elements of said vector MX₀ where M is said oversampling factor employed in said FFT calculation. The above mentioned Steps 2.3.1 - 2.3.2 of Fig. 7 can be replaced by an alternative method including (see Steps 2.4.1 - 2.4.2 of Fig. 8):

- (Step 2.4.1) said iterative procedure includes a preliminary initialization comprising: a) setting the iteration index i to 1, the initial estimate of A equal to said vector element (X_0jg)_/ the initial estimate of A = 27T5 F_DFT to zero and the initial estimate of F to said main portion (a F_DFT ) ; b) computing the real quantity p = F/F_DFT; then iteratively

- (Step 2.4.2) the i-th iteration is fed by the estimates of D, F and A and the real quantity p computed in the previous iteration, to generate new estimates of D, F, A and p, for a predetermined number of iterations, wherein each iteration includes: a) a first step where the residual D is calculated solving the first order equation b^A + % = 0 OR the second order equation a_dD² F ί¾D F CQ ⁼ 0 in the variable D with: where ¾{x} denotes the real part, 3{x} denotes the imaginary part of the complex quantity x and {X_kp; k = 1,2,3} is computed as: where x_cn is the n -th element of the said complex; said coefficients ({X_kp; k = 1,2,3}) can be also computed by interpolating the elements of said spectrum (X_fc); discarding one of the two solutions, which does not satisfy the inequalities: ifa = 0 otherwise if said second order equation is solved; then the new estimate of F is obtained by adding said residual (A F_dft) to said estimate of F; b) a second step wherein the new estimate of p is computed by dividing said new estimate of F by F_dft; c) a third step wherein said new estimate of F is employed to compute the new estimate of A as A = X_int generated by interpolating the elements of said vector MX₀ where M is said oversampling factor employed in said FFT calculations. The above-described methods for the estimation of the frequency, amplitude and phase of a complex sequence can be employed for range (R) calculation of a point target detected through a single-input single-output frequency modulated continuous wave radar system as where m is the radar chirp rate, c is the speed of light, f_s is the sampling frequency and F is said estimated normalized frequency calculated according to said methods.

Additionally, said methods can be employed for angle of arrival (Q) calculation of a point target having range (R) detected through a multiple-input multiple-output frequency modulated continuous wave radar system equipped with a plurality of transmitting (TX) and receiving (RX) antennas, wherein each couple of said TX and RX antennas is replaced with the equivalent virtual antenna; said virtual antennas forming a virtual array composed by one uniform linear array; said angle of arrival being computed as: q= arcsin here, l is the radar wavelength, d is the distance between two adjacent virtual antennas of said uniform linear array, and F is said estimated normalized frequency according to said methods. Estimation of single delays

All the estimation techniques developed in the previous two paragraphs refer to the real (complex) sequence {x_r,n} ({x_c,n})_r whose n-th element is expressed by eq. (10) (eq. (15)). As shown in the previous section, these models are suitable to represent the sampled downconverted signal made available by the receiver of a colocated FMCW MIMO radar. However, similar estimation methods can be also developed for a colocated SFCW MIMO radar system, since the complex signal model (18) illustrated for this system is very similar to the model (15) referring to its complex FMCW counterpart, the only differences being represented by the fact that: a) the physical meaning of the parameter F_; , since this represents a normalized delay in eq. (18) and a normalized frequency in eq. (15); b) the sign of the argument of all the complex exponentials appearing in the RHS of eq. (18), since this is the opposite of that of the corresponding functions contained in the RHS of eq. (15). The implications of this sign reversal can be summarised as follows. In the derivation of the single tone estimator, the function X(F) (38) is replaced by where x_cn is the n-th sample of the received signal sequence (see eq. (18)). Consequently, the vectors X₀ (42) and X_¾ (70) are generated by taking the iV₀ order inverse DFT (IDFT) of the zero- padded vectors x_{0 ZP} (43) and x_{fc ZP} (71), respectively; note that the n-th element of the N₀-th dimensional vector x_0,zp (^x _/c_,zp ) is equal to x_cn (n^{k ■} x_cn; see eq. (69)) for n = 0 , 1, ..., N — 1 and is equal to zero for n > N . Two algorithms for estimating a single delay can be easily developed by following the same line of reasoning as the CSFEs; these are called complex single delay estimators (CSDEs). More specifically, the CSDEs compute an estimate d of the delay residual d by solving eq. (74) in the variable D (57); in this case, the coefficients a(p) and c(p) are still expressed by eqs. (105) and (107), respectively, whereas b(p)A -¾{i^*X_2,p}, (109) where A is the estimate of the complex gain A computed on the basis of eq. (102). Moreover, in the evaluation of the coefficients {a(p), b(p), c(p)} of eq. (74), the quantity X_kp(i-q (see eq. (68)) is replaced by where pO^-1-¹ is still defined by eq. (94).

Preferably, the method for the estimation of the delay, the amplitude and the phase of a complex sequence (see Steps 2.5.0

- 2.5.2 of Fig. 9) can be described as follows:

- (Step 2.5.0) representation of a tone normalized delay T = t Af , where Af is the frequency step size of the employed radar, as sum of a main portion (QCF_idft) and a residual portion (5F_idft), where F_IDFT = l/iV₀, where iV₀ is the IFFT order, wherein said main portion is determined by searching the index (a) of the vector element X_oa corresponding to the maximum absolute value of the same vector X₀ and said residual portion is found through an iterative procedure together with the corresponding target complex amplitude A = a exp(-jxp), where a is said target amplitude and y is said target phase.

- (Step 2.5.1) said iterative procedure includes a preliminary initialization comprising: setting the iteration index i to 1, the initial estimate of A equal to said vector element (Xo_,a) _r the initial estimate of D = 2i F_IDFT to zero and the initial estimate of T to said main portion (aF_IDFT); then iteratively

- (Step 2.5.2) the i-th iteration is fed by the estimates of D, T and A computed in the previous iteration, to generate new estimates of D , T and A , for a predetermined number of iterations, wherein each iteration includes a) a first step where the residual D is calculated solving the first order equation Z¾D + c¾ = 0 OR the second order equation a_dD² + Z¾D + CQ = 0 in the variable D with: where ¾{x} denotes the real part and 3{x} denotes the imaginary part of the complex quantity x ; discarding one of the two solutions, which does not satisfy the inequalities: ifa = 0 otherwise if said second order equation is solved; then the new estimate of T is obtained by adding said residual (A F_idft) to said main portion (af_IDFT); b) a second step wherein said new estimate of T is employed to compute the new estimate of A as A = X_int generated by interpolating the elements of said vector MX₀ where M is said oversampling factor employed in said IFFT calculation. The above mentioned Steps 2.5.1 - 2.5.2 of Fig. 9 can be replaced by an alternative method including (see Steps 2.6.1 - 2.6.2 of Fig. 10):

- (Step 2.6.1) said iterative procedure includes a preliminary initialization comprising: a) setting the iteration index i to 1, the initial estimate of A equal to said vector element X_0¾, the initial estimate of D = 2ir<5F_IDFT to zero and the initial estimate of T to said main portion (af_IDFT); b) computing the real quantity p = T/F_IDFT; then iteratively - (Step 2.6.2) the i-th iteration is fed by the estimates of D,

T, A and the real quantity p computed in the previous iteration, to generate new estimates of D, T, A and p, for a predetermined number of iterations, wherein each iteration includes: a) a first step where the residual D is calculated solving the first order equation b_%A + c¾ = 0 OR the second order equation ¾D² + ¾D + ¾ = 0 in the variable D with: where ¾{x} denotes the real part, 3{x} denotes the imaginary part of the complex quant computed as: where x_{c n} is the n-th element of the said complex signal; said coefficients ({X_kp; k = 1,2,3}) can be also computed by interpolating the samples of said impulse response (X_¾); discarding one of the two solutions, which does not satisfy the inequalities: if a = 0 otherwise if said second order equation is solved; then the new estimate of T is obtained by adding said residual (AF_idft) to said estimate of T ; b) a second step wherein the new estimate of p is computed by dividing said new estimate of T by F_idft; c) a third step wherein said new estimate of T is employed to compute the new estimate of is generated by interpolating the elements of said vector MX₀ where M is said oversampling factor employed in said IFFT calculations.

The above-described methods for the estimation of the delay, amplitude and phase of a complex sequence can be employed for range (R) calculation of a point target detected through a single-input single-output stepped frequency continuous wave radar system as where D/ is the frequency step size of the employed radar system, c is the speed of light and T is said estimated normalized delay, calculated according to said methods.

Additionally, said methods can be employed for angle of arrival (Q) calculation of a point target having range (R) detected through a multiple-input multiple-output stepped frequency continuous wave radar system equipped with a plurality of transmitting (TX) and receiving (RX) antennas, wherein each couple of said TX and RX antennas is replaced with the equivalent virtual antenna; said virtual antennas forming a virtual array composed by one uniform linear array; said angle of arrival being computed as: here, l is the radar wavelength, d is the distance between two adjacent virtual antennas of said uniform linear array, and T is said estimated normalized delay according to said methods.

This invention can be implemented advantageously in a computer program comprising program code means for performing one or more steps of such method, when such program is run on a computer. For this reason, the patent shall also cover such computer program and the computer-readable medium that comprises a recorded message, such computer-readable medium comprising the program code means for performing one or more steps of such method, when such program is run on a computer, such as a computing device associated to a radar system or to a base station of cellular telecommunication, such as 5G and similar ones.

Many changes, modifications, variations and other uses and applications of the subject invention will become apparent to those skilled in the art after considering the specification and the accompanying drawings which disclose preferred embodiments thereof as described in the appended claims.

The features disclosed in the prior art background are introduced only in order to better understand the invention and not as a declaration about the existence of known prior art. In addition, said features define the context of the present invention, thus such features shall be considered in common with the detailed description.

Further implementation details will not be described, as the man skilled in the art is able to carry out the invention starting from the teaching of the above description.

Claims

1. Method for computer estimation of the frequency, amplitude and phase of a single tone defining a signal acquired by means of a receiver, wherein the tone is real or complex; said method including:

- (Step 1: FFT processing) estimation of a spectrum (X₀) of said signal and its first N_FFT — 1 derivatives ({X_fc;k= 1,.. N_{FFT —} 1}) with N_FFT = 3 when said signal is real and N_FFT = 4 when said signal is complex, through FFT calculation with an oversampling factor (M),

- (Step 2) iterative calculation of an estimate of the parameters of the most dominant tone, including phase (y), amplitude (a) and frequency (/) on the basis of said spectrum, and its derivatives, through N_{FFT —} 1 FFTs.

2. [SFE] Method according to claim 1, wherein said signal is real, said (Step 2) calculation includes:

- (Step 2.1.0) representation of a tone normalized frequency F = f/f_sr where f_s is the sampling frequency, as sum of a main portion (aF_dft) and a residual portion (5F_dft), where F_DFT = l/iV₀, where iV₀ is the FFT order, wherein said main portion is determined by searching the index (a) of the vector element X_oa corresponding to the maximum absolute value in a first half plus one elements of the same vector X₀ and said residual frequency is found through said iterative calculation together with the corresponding tone complex amplitude C = aexp(jxp)/2, where a is said tone amplitude and y is said tone phase.

3. [SFE#1] Method according to claim 2, wherein said iterative procedure includes a preliminary initialization (Step 2.1.1) comprising: a) setting the iteration index i to 1, the initial estimate of tone complex amplitude C equal to said vector element (X_0,a)_t the initial estimate of A = 2pd F_DFT to zero and the initial estimate of F to said main portion (CCF_dft); b) computing the complex coefficients (K_p(2a); p = 1,2,3} as the 2a- th element of the iV₀-order FFT of the sequence n^p, with p = 1, 2, 3 ; then iteratively

- (Step 2.1.2) the i-th iteration is fed by the estimates of D, F and C computed in the previous iteration, to generate new estimates of D , F and C , for a predetermined number of iterations, wherein each iteration includes a) a first step where the residual D is calculated solving the first order equation b^K + c¾ = 0 OR the second order equation a_dD² + ί¾D + C_Q ⁼ 0 in the variable D with: where ¾{x} denotes the real part and 3{x} denotes the imaginary part of the complex quantity x; discarding one of the two solutions, which does not satisfy the inequalities: if a = 0 otherwise if said second order equation is solved; then the new estimate of F is obtained by adding said residual (AF_dft) to said main portion (aF_DFT); b) a second step wherein said new estimate of F is employed to compute the new estimate of C as where is generated by interpolating the elements of said vector MX₀ where M is said oversampling factor employed in said FFT calculations.

4. [SFE#2] Method according to claim 2, wherein said iterative procedure includes a preliminary initialization (Step 2.2.1) comprising: a) setting the iteration index i to 1, the initial estimate of C equal to said vector element (X_0,a)_t the initial estimate of A =

2i F_DFT to zero and the initial estimate of F to said main portion (af_DFT); b) computing the complex coefficients (K_p(2a);p = 1,2,3} as the 2a- th element of the iV₀-order FFT of the sequence n^p , with p = 1,2,3; c) computing the real quantity p=F/F_DFT; then iteratively

- (Step 2.2.2) the i-th iteration is fed by the estimates of D, F and C and the real quantity p computed in the previous iteration, to generate new estimates of D, F, C and p, for a predetermined number of iterations, wherein each iteration includes a) a first step where the residual D is calculated solving the first order equation Z¾D + c¾ = 0 OR the second order equation a_dD² + Z¾D + C_Q ⁼ 0 in the variable D with: a_? - 23 {(r)² K₃(2p)j, where ¾{x} denotes the real part, 3{x} denotes the imaginary part of the complex quantity x, {K_r(2b); p = 1,2,3} are computed as: with p = 1, 2, 3 as: where c_Lh is the n-th element of the said real signal; said coefficients ({X_k,p> k = 1,2}) are also computed by interpolating the elements of said spectrum (X_fc ); discarding one of the two solutions, which does not satisfy the inequalities: if a = 0 otherwise if said second order equation is solved; then the new estimate of F is obtained by adding said residual (A F_dft) to said estimate of F; b) a second step wherein the new estimate of p is computed by dividing said new estimate of F by F_dft; c) a third step wherein said new estimate of F is employed to compute the new estimate of C as: where is generated by interpolating the elements of said vector MX₀ where M is said oversampling factor employed in said FFT calculations.

5. [CSFE] Method according to claim 1, wherein said signal is complex, said (Step 2) calculation includes:

- (Step 2.3.0) representation of a tone normalized frequency F = f/f_sr where f_s is the sampling frequency, as sum of a main portion (aF_dft) and a residual portion (5F_dft), where F_DFT = l/iV₀, where N₀ is the FFT order, wherein said main portion is determined by searching the index (a) of the vector element X_oa corresponding to the maximum absolute value of the same vector X₀ and said residual frequency is found through said iterative calculation together with the corresponding tone complex amplitude A = aexp(jxp), where a is said target amplitude and y is said target phase.

6. [CSFE#1] Method according to claim 5, wherein said iterative procedure includes a preliminary initialization (Step 2.3.1) comprising: setting the iteration index i to 1, the initial estimate of A equal to said vector element (Xo_,a)t the initial estimate of A = 2ir<5F_DFT to zero and the initial estimate of F to said main portion (af_DFT); then iteratively

- (Step 2.3.2) the i-th iteration is fed by the estimates of D, F and A computed in the previous iteration, to generate new estimates of D , F and A , for a predetermined number of iterations, wherein each iteration includes a) a first step where the residual D is calculated solving the first order equation b^L·+ c¾ = 0 OR the second order equation a_dD² + Z¾D + C_Q ⁼ 0 in the variable D with: where ¾{x} denotes the real part and 3{x} denotes the imaginary part of the complex quantity x; discarding one of the two solutions, which does not satisfy the inequalities: ifa = 0 otherwise if said second order equation is solved; then the new estimate of F is obtained by adding said residual (AF_dft) to said main portion (af_DFT); b) a second step wherein said new estimate of F is employed to compute the new estimate of A as A = generated by interpolating the elements of said vector MX₀ where M is said oversampling factor employed in said FFT calculation.

7. [CSFE#2] Method according to claim 5, wherein said iterative procedure includes a preliminary initialization (Step 2.4.1) comprising: a) setting the iteration index i to 1, the initial estimate of A equal to said vector element (X_0@) , the initial estimate of A = 2i F_DFT to zero and the initial estimate of F to said main portion (aF_DFT); b) computing the real quantity p=F/F_DFT; then iteratively

- (Step 2.4.2) the i-th iteration is fed by the estimates of D, F and A and the real quantity p computed in the previous iteration, to generate new estimates of D, F, A and p, for a predetermined number of iterations, wherein each iteration includes: a) a first step where the residual D is calculated solving the first order equation b^A+ C_Q = 0 OR the second order equation where ¾{x} denotes the real part, denotes the imaginary part of the complex quantity x and {X_kp;k = 1,2,3} is computed as: where x_cn is the n -th element of the said complex; said coefficients ({X_kp; k = 1,2,3}) can be also computed by interpolating the elements of said spectrum (X_fc); discarding one of the two solutions, which does not satisfy the inequalities: ifa = 0 otherwise if said second order equation is solved; then the new estimate of F is obtained by adding said residual (A F_dft) to said estimate of F; b) a second step wherein the new estimate of p is computed by dividing said new estimate of F by F_dft; c) a third step wherein said new estimate of F is employed to compute the new estimate of A as A = X_int is generated by interpolating the elements of said vector MX₀ where M is said oversampling factor employed in said FFT calculations.

8. Method for range (R) calculation of a point target detected through a single-input single-output frequency modulated continuous wave radar system as where m is the radar chirp rate, c is the speed of light, f_s is the sampling frequency and F is said estimated normalized frequency calculated according to any one of previous claims from 1 to 7.

9. Method for angle of arrival (Q) calculation of a point target having range (R) detected through a multiple-input multiple- output frequency modulated continuous wave radar system equipped with a plurality of transmitting (TX) and receiving (RX) antennas, wherein each couple of said TX and RX antennas is replaced with the equivalent virtual antenna; said virtual antennas forming a virtual array composed by one uniform linear array; said angle of arrival being computed as:

Q= arcsin here, l is the radar wavelength, d is the distance between two adjacent virtual antennas of said uniform linear array, and F is said estimated normalized frequency according to any one of the claims from 1 to 7.

10. Method for computer estimation of the delay, amplitude and phase of a single complex tone defining a signal acquired by means of a receiver; the method including:

- (Step 1: IFFT processing) estimation of the impulse response (X₀) characterizing the communication channel and its first

— (v) N_FFT — 1 derivatives ({X_fc ;k = 1,.. N_FFT — 1}) with N_FFT = 4, through IFFT calculation, with an oversampling factor (M ),

- (Step 2) iterative calculation of an estimate of the parameters of the most dominant tone, including phase (y), amplitude (a) and delay (t) on the basis of said impulse response, and its derivatives, through N_FFT — 1 IFFTs.

11. Method according to claim 10, wherein said (Step 2) calculation includes: - (Step 2.5.0) representation of a tone normalized delay T = t Af , where Af is the frequency step size of the employed radar, as sum of a main portion (OCF_idft) and a residual portion (5F_idft), where F_IDFT = l/iV₀, where iV₀ is the IFFT order, wherein said main portion is determined by searching the index (a) of the vector element X_oa corresponding to the maximum absolute value of the same vector X₀ and said residual portion is found through said iterative calculation together with the corresponding target complex amplitude A = a exp(-jxp), where a is said target amplitude and y is said target phase.

12. [CSDE#1] Method according to claim 11, wherein said iterative procedure includes a preliminary initialization (Step 2.5.1) comprising: setting the iteration index i to 1, the initial estimate of A equal to said vector element (Xo_,a) _/ the initial estimate of A = 27Z1>F_IDFT to zero and the initial estimate of T to said main portion (aF_IDFT); then iteratively

- (Step 2.5.2) the i-th iteration is fed by the estimates of D, T and A computed in the previous iteration, to generate new estimates of D , T and A , for a predetermined number of iterations, wherein each iteration includes a) a first step where the residual D is calculated solving the first order equation OR the second order equation in the variable D with: where ¾{x} denotes the real part and 3{x} denotes the imaginary part of the complex quantity x ; discarding one of the two solutions, which does not satisfy the inequalities: if a = 0 otherwise if said second order equation is solved; then the new estimate of T is obtained by adding said residual to said main portion (aF_IDFT); b) a second step wherein said new estimate of T is employed to compute the new estimate of and is generated by interpolating the elements of said vector where M is said oversampling factor employed in said IFFT calculation.

13. [CSDE#2] Method according to claim 11, wherein said iterative procedure includes a preliminary initialization (Step 2.6.1) comprising: a) setting the iteration index i to 1, the initial estimate of A equal to said vector element the initial estimate of D = 2i F_IDFT to zero and the initial estimate of T to said main portion b) computing the real quantity then iteratively

- (Step 2.6.2) the i-th iteration is fed by the estimates of D, T, A and the real quantity p computed in the previous iteration, to generate new estimates of D, T, A and p, for a predetermined number of iterations, wherein each iteration includes: a) a first step where the residual D is calculated solving the first order equation OR the second order equation : where ¾{x} denotes the real part, 3{x} denotes the imaginary part of the complex quantity x and {X_kp;k = 1,2,3} is computed as: where x_cn is the n-th element of the said complex signal; said coefficients ({X_k,p>k = 1,2,3}) can be also computed by interpolating the samples of said impulse response (X_fc); discarding one of the two solutions, which does not satisfy the inequalities: ifa = 0 otherwise if said second order equation is solved; then the new estimate of T is obtained by adding said residual to said estimate of T; b) a second step wherein the new estimate of p is computed by dividing said new estimate of T by F_idft; c) a third step wherein said new estimate of T is employed to compute the new estimate of A as A = X_int generated by interpolating the elements of said vector where M is said oversampling factor employed in said IFFT calculations.

14. Method for range calculation of a point target detected through a single-input single-output stepped frequency continuous wave radar system as where D/ is the frequency step size of the employed radar system, c is the speed of light and T is said estimated normalized delay, calculated according to any one of the previous claims from 10 to 13.

15. Method for angle of arrival (Q) calculation of a point target having range ( R) detected through a multiple-input multiple- output stepped frequency continuous wave radar system equipped with a plurality of transmitting (TX) and receiving (RX) antennas, wherein each couple of said TX and RX antennas is replaced with the equivalent virtual antenna; said virtual antennas forming a virtual array composed by one uniform linear array; said angle of arrival being computed as: here, l is the radar wavelength, d is the distance between two adjacent virtual antennas of said uniform linear array, and T is said estimated normalized delay according to any one of the claims 10 - 13.

16. A MIMO FMCW radar system including:

- a MIMO FMCW radar equipped with a plurality of transmitting (TX) and receiving (RX) antennas, arranged to generate real or complex signals in response to a propagation scenario including a plurality of point targets, and

- computing means arranged to receive signals generated by said transmitting and receiving antennas and configured to replace each couple of said transmitting and receiving antennas with an equivalent model defining virtual antennas and to perform all the steps of any one of the previous claims from 1 to 9.

17. A SISO FMCW radar system including:

- a SISO FMCW radar equipped with a single transmitting (TX) and receiving (RX) antennas, arranged to generate real or complex signals in response to a propagation scenario including a plurality of point targets, and

- computing means arranged to receive signals generated by said receiving antenna and configured to perform all the steps of any one of the previous claims from 1 to 9.

18. A MIMO SFCW radar system including:

- a MIMO SFCW radar equipped with a plurality of transmitting (TX) and receiving (RX) antennas, arranged to generate complex signals in response to a propagation scenario including a plurality of point targets, and

- computing means arranged to receive signals generated by said transmitting and receiving antennas and configured to replace each couple of said transmitting and receiving antennas with an equivalent model defining virtual antennas and to perform all the steps of any one of the previous claims from 10 to 15.

19. A SISO SFCW radar system including:

- a SISO SFCW radar equipped with a single transmitting (TX) and receiving (RX) antennas, arranged to generate complex signals in response to a propagation scenario including a plurality of point targets, and - computing means arranged to receive signals generated by said receiving antenna and configured to perform all the steps of any one of the previous claims from 10 to 15.