Millimeter-wave Mobile Sensing and Environment Mapping: Models, Algorithms and Validation

Integrating efficient connectivity, positioning and sensing functionalities into 5G New Radio (NR) and beyond mobile cellular systems is one timely research paradigm, especially at mm-wave and sub-THz bands. In this article, we address the radio-based sensing and environment mapping prospects with specific emphasis on the user equipment (UE) side. We first describe an efficient l1-regularized least-squares (LS) approach to obtain sparse range-angle charts at individual measurement or sensing locations. For the subsequent environment mapping, we then introduce a novel state model for mapping diffuse and specular scattering, which allows efficient tracking of individual scatterers over time using interacting multiple model (IMM) extended Kalman filter and smoother. Also the related measurement selection and data association problems are addressed. We provide extensive numerical indoor mapping results at the 28 GHz band deploying OFDM-based 5G NR uplink waveform with 400 MHz channel bandwidth, covering both accurate ray-tracing based as well as actual RF measurement results. The results illustrate the superiority of the dynamic tracking-based solutions, compared to static reference methods, while overall demonstrate the excellent prospects of radio-based mobile environment sensing and mapping in future mm-wave networks.

sensing capabilities to the networks and the corresponding user devices [6]- [10]. To this end, the development of joint communication and sensing (JCAS) systems that can perform both functionalities while sharing the same transmit waveforms and hardware platforms is a very timely research area [11]- [14], with applications, e.g., in vehicular systems and indoor mapping [11], [15]- [18]. The large channel bandwidths available at mmwave frequencies, together with highly directional antenna arrays, form the technical basis for high-accuracy radio positioning and sensing [4], [19], [20]. While 3GPP 5G NR supports currently channel bandwidths up to 400 MHz [21], further increased channel bandwidths are expected through the NR evolution towards the sub-THz regime [22], paving eventually the way to future 6G networks [16], [23].
In general, in typical outdoor mm-wave scenarios with densely deployed network nodes or radio heads, the user equipment (UE) is likely to have line-of-sight (LoS) connection to the network while also being within the coverage range of multiple network nodes -both being issues that further facilitate high-accuracy positioning [2]- [4]. However, in indoor deployment scenarios, the propagation characteristics and thus the positioning become much more challenging due to the high path loss as well as the reduced penetration and diffraction [24], [25]. To address such challenges, the concept of simultaneous localization and mapping (SLAM) has been studied. To this end, in [26]- [36], different SLAM algorithms aim to estimate the UE position and the environment scatterer positions by exploiting the multipath component (MPC) information from physical or fixed anchors, e.g., 5G NR base stations.
In [34], a map-free indoor localization method using ultrawideband (UWB) mm-wave signals is proposed, utilizing the estimated scatterer positions of the environment as virtual anchors to obtain the UE location. However, a snap-shot approach is taken, thus ignoring the correlation in sequences of measurements obtained by moving sensors. The work in [35], in turn, described a SLAM architecture for mm-wave localization through obstacle detection and dimensioning for indoor environments. Furthermore, a message passing-based estimator which jointly estimates the position and orientation of the UE while sensing the environment's reflectors is presented in [36]. The works [35] and [36] do not consider the actual MPC estimation and do not take into account the different nature of the specular reflections and diffuse scattering in the actual mapping task.
In contrast to the traditional SLAM methods which fuse information from fixed anchors and moving agents, alternative techniques consider user-centric approaches where localization and mapping functionalities are implemented at each agent independently, without the need for fixed infrastructure [37]. In this case, each moving agent is equipped with transmit and receive antennas, similar to a monostatic radar system, that measures the MPCs and performs SLAM accordingly.
In [37], it is assumed that a reference map or floorplan is available, hence focusing mostly on UE localization, while also a recursive update method to account for the uncertainties of the floorplan in user-centric SLAM system is proposed. Only ranging or delay measurements are considered in [37] while the environment related state model comprises essentially of wall segments. In [19], a personal mobile radar concept for environment sensing or mapping is studied, building on large antenna arrays at mm-waves. The mapping related state variables comprise a collection of radar cross-section (RCS) values for a predefined grid of geographical cells covering uniformly the overall considered area. Such an approach leads easily to thousands of state variables, and hence large computational complexity. In [38], a JCAS system is pursued for smart home scenarios, investigating different integration approaches for traditional sensing and WiFi devices. In general, due to the notable challenges in performing RF measurements at mm-wave bands, only very few works have supported their positioning/sensing/mapping solutions with empirical results to assess and demonstrate the performance with real mm-wave transceiver and antenna hardware [25], [34], [39]- [41].
In this article, we describe a novel application framework for user-centric mm-wave indoor mapping systems, propose associated signal processing methods in 5G NR UE context, assuming a known agent trajectory, and perform validation with ray-tracing and RF measurements. In the considered approach, illustrated conceptually in Fig. 1, the UE senses the surrounding environment through orthogonal frequency-division multiplexing (OFDM)-based beamformed uplink transmissions while then observing and collecting the target reflections with synchronous receiver beamforming patterns. This is followed by range-angle processing, while the corresponding rangeangle charts are then further post-processed in the subsequent dynamic mapping stage. The considered radio-based sensing and mapping concept exploits the UE mobility, which collects correlated sensing measurements at different locations to better reconstruct the environment. We specifically address the moving nature of the sensing and measurement device, in the form of advanced tracking-based dynamic mapping solutions, incorporating both specular reflections and diffuse scattering. With the proposed tracking approach, besides mapping of fixed structures, it is in general possible to track moving targets, adapt to the changes in the environment, as well as predict the mobility of a tracked target.
With respect to the overall state-of-the-art literature in the field, the contributions and novelty of the article can be summarized as follows: • An efficient 1 -regularized least-squares (LS) formulation to obtain sparse range-angle charts at individual measurement or sensing locations, expressed for OFDM systems, is provided; • A novel state model for mapping diffuse and specular scattering is introduced, which allows efficient tracking of individual scatterers over time using interacting multiple model (IMM) extended Kalman filter/smoother; • Multi-bounce or double-reflection challenge of practical complex environments is addressed through the available received signal strength (RSS) measurements and applicable pathloss models, and corresponding novel measurement selection and data association methods are devised; • Comprehensive numerical indoor mapping results at 28 GHz band deploying OFDM-based 5G NR uplink waveform with 400 MHz channel bandwidth are provided, covering both accurate ray-tracing based as well as actual RF measurement results. The rest of this article is organized as follows. The fundamental system model for collecting beamformed sensing measurements and the associated measurement system geometry are described in Section II. The LS-based rangeangle processing methods for individual sensing locations are addressed and derived in Section III, while Section IV presents the proposed tracking-based dynamic mapping solutions for a sequence of sensing locations, respectively. In Section V, the evaluation scenario, ray-tracing environment and the experimental RF measurement environment and equipment are described. Obtained numerical results are provided and analyzed in Section VI, while the conclusions are provided in Section VII. Finally, further details of the LS-based rangeangle processing solution are described in the Appendices A and B.

A. Basics
The sensing and mapping functionality of the NR UE builds on the OFDM-based uplink transmit waveform whose complex I/Q samples are fully known in the device. This is conceptually illustrated in Fig. 1(a). Furthermore, the mm-wave UE is assumed to be equipped with a directive antenna system for both transmit (TX) and receive (RX) functions, which allows for accurate angular measurements to map the environment. In the following system model description, for notational convenience, uniform linear arrays (ULAs) are considered for both the TX and RX sides. However, it is noted for clarity that directive horn antennas are used in the later RF measurements due to the available hardware limitations. Additionally, in this work, we specifically focus on sensing the azimuth domain to create a 2D map of the environment. However, the same principles can be applied for sensing in the elevation domain and subsequently for extracting 3D maps.
For given sensing direction, we assume that the UE transmits a sequence of M beamformed OFDM symbols on N active subcarriers, with x n,m denoting the transmitted frequencydomain symbol at n th subcarrier in m th OFDM symbol. The presented system model assumes a static channel across the M transmitted symbols which can be extended to consider moving targets as shown in [13], [42]. To this end, the frequencydomain OFDM-based radar model with K targets or MPCs reads [13], [42], [43] where v n,m , w TX ∈ C NTX×1 and w RX ∈ C NRX×1 are the complex Gaussian noise and the TX and RX array beamforming weights, respectively. Moreover, N TX and N RX denote the number of antenna elements in the TX and RX array, and A TX,n ∈ C NTX×K and A RX,n ∈ C NRX×K are the steering matrices for the TX and RX array, respectively. The multipath coefficients of a sensed environment with K target reflections are modelled through Γ n = diag(γ 0,n , . . . , γ K−1,n ) ∈ C K×K that is a diagonal matrix with its k th element being of the form In above, τ k , d k , σ RCS k and f c model the k th target's propagation delay, propagation distance, RCS and the system carrier frequency, respectively, stemming from the well-known radar range equation [42]. furthermore, λ n refers to the wavelength of the n th subcarrier while ∆f denotes the subcarrier spacing of the OFDM waveform.
Under the ULA assumption, the TX array steering vector for a k th target with azimuth angle ϕ TX k can be expressed at subcarrier n as where Υ n (ϕ TX k ) = 2π υ λn sin (ϕ TX k ) is the corresponding electrical angle of departure (AoD) for the k th target, wherein υ is the antenna element separation. The RX steering vector a RX,n (ϕ RX k ) and the corresponding electrical angle of arrival (AoA) Υ n (ϕ RX k ) are obtained similarly. As a result, the set of steering vectors for the considered K target reflections can be expressed with compact matrix notation as A TX,n = [a TX,n (ϕ RX 0 ), . . . , a TX,n (ϕ RX K−1 )] and A RX,n = [a RX,n (ϕ RX 0 ), . . . , a RX,n (ϕ RX K−1 )] which are deployed in (1).

B. Geometry
Next, we consider the specific characteristics related to the fundamental system geometry shown in Fig. 2 where the sensing device is located at position p UE l = (x UE l , y UE l ) with l denoting the sensing instance or location index. Later in this work, in Section IV, we consider the fact that the UE is moving, hence leading to a sequence of sensing locations p UE l = (x UE l , y UE l ), l = 1, ..., L. However, here and in the following Section III, an arbitrary given sensing location is considered.
Like illustrated in the figure, we allow for separate TX and RX arrays in our system model, to relax the self-interference challenge in OFDM radars [13], [44]. Therefore, with separate TX array position p TX l = (x TX l , y TX l ) and RX array position p RX l = (x RX l , y RX l ), the physical propagation delay of the k th target reflection τ l,k = (d TX l,k + d RX l,k )/c 0 is determined cumulatively according to the distance between the TX array and target position d TX l,k , and the distance between the RX array and target position d RX l,k . Consequently, the distance (i.e., range) between the UE position p UE l and the k th target position p SC l,k with azimuth angle ϕ UE l,k can be expressed as where c 0 is the speed of light and d ant corresponds to the separation between TX and RX arrays. For the special case of co-located or joint TX-RX antenna system, d ant = 0 while the above expression is written in a general form. Furthermore, for a given physical propagation delay τ l,k and the system geometry shown in Fig. 2, the transmit angle ϕ TX l,k -and similarly the receive angle -can be described as where ρ(τ l,k , ϕ UE l,k ) is defined in (4) and atan2(y, x) is a fourquadrant inverse tangent function.
Assuming then that the distances of the considered scatterers satisfy the far-field condition, we can express the TX and RX gain patterns associated to a target or scatterer k as For simplicity, in the rest of the paper, we assume frequency-flat gain patterns with g TX,n (ϕ TX l,k ) = g TX (ϕ TX l,k ) and g RX,n (ϕ RX l,k ) = g RX (ϕ RX l,k ). Finally, the combined TX-RX antenna pattern is defined as For generality, we also note that in the JCAS literature, different approaches have been considered to optimize the TX and RX beamforming weights, w TX and w RX , see, e.g., [45]- [48]. Specifically, different advanced techniques enabling the transmission of multiple simultaneous beams have been proposed in these recent works. In this paper, however, we consider more ordinary single-beam approach due to its implementation simplicity, while defer the potential multi-beam considerations to our follow-up future work.

III. RANGE-ANGLE PROCESSING
In this section, we start by extending the basic signal model in (1) to the case of multiple sensing directions for rangeangle charting. Then, we introduce a novel formulation of the OFDM radar range-angle processing problem for the proposed environment mapping scenario and propose an 1 -regularized LS estimation problem to obtain sparse range-angle charts to efficiently facilitate subsequent mapping phases, while also shortly quantify the involved processing complexity. It is noted that the presentation in this section considers an individual arbitrary sensing location p UE l . For readers' convenience, the main variables used in this section are summarized in Table I.  where denotes the Hadamard (element-wise) product, C R and C ϕ are the numbers of range and azimuth cells in the discretized range-angle chart, respectively, and b p,q is the reflection coefficient at the (p, q) th range-angle cell located at (τ p ,φ q ). The variables g(ϕ)

A. Problem Statement
and refer to the combined TX-RX antenna pattern, defined in (7), and the frequency-domain additive noise vector, respectively. Furthermore, in (8), the frequency-domain steering vector for a specific cell's delayτ p reads The purpose of the range-angle processing is to estimate the reflection coefficients b p,q = (B) p,q at predefined rangeangle grid locations (τ p ,φ q ) with p = 0, . . . , C R − 1 and q = 0, . . . , C ϕ −1, while exploiting the knowledge of the combined antenna pattern. Note that b p,q = 0 if there is no scatterer at the (p, q) th range-azimuth cell. For notational convenience, we can next rewrite the frequency-domain received signal or measurement vector (8) as Specifically, the matrix B defined in (11d) represents the rangeangle chart to be estimated. Aggregating (10) over all I sensing directions, the observation matrix for the m th OFDM symbol at a given UE location can be expressed as At any given l th UE location, the range-angle estimation problem can then be formally defined as follows. Given the transmit symbols {X m } M −1 m=0 , the steering matrix C (known through the delay grid locations {τ p } CR−1 p=0 ), the antenna pattern matrix G (known through the angular grid locations {φ q } Cϕ−1 q=0 and the sensing directions (12).

B. Regularized LS for Sparse Range-Angle Processing
To harness the sparsity of the mm-wave propagation channels, and thus to obtain sparse range-angle charts from (12), we propose to consider the following 1 -regularized LS problem: where , Φ is a dictionary matrix defined in Appendix A, and λ denotes the regularization parameter. Here, vec (·) denotes the vectorization operator. For reference, Appendix A also describes the standard LS solution B LS of (14) without regularization. To solve the sparse map recovery problem in (14), we resort to the iterative shrinkage/thresholding algorithm (ISTA) [49]- [51], as outlined in Algorithm 1, where λ max (·) denotes the maximum eigenvalue of a matrix and (·) + yields the positive part of a real number while sgn(z) = z/|z| yields the sign of a complex number.
The obtained sparse range-angle charts are then subject to a target detection phase that provides the input for the subsequent mapping phases, described in Section IV. Specifically, a local maximum detection is implemented to the ISTA output, considering a minimum target separation in both range and angle domains. In addition, we limit the maximum number of detected targets per sensing location to maintain the complexity of the following tracking-based dynamic mapping phase. Moreover, a threshold test is considered to discard weak and thereon irrelevant target reflections. Further specifics in terms of the numerical threshold values are provided along the numerical results in Sections V and VI.
Finally, the corresponding Cartesian coordinates can be calculated for any given point of the range-angle chart as follows. Considering a UE location p UE Fig. 2, the Cartesian coordinates of the (p, q) th range-angle cell located at (τ p ,φ q ) can be calculated as where a = {p + qC R + lC R C ϕ , ∀p, ∀q, ∀l}, and ρ(τ, ϕ) is defined in (4). The parameter ϕ OR l denotes the UE orientation angle as illustrated in Fig. 2.
until convergence

C. Complexity Analysis
The computational complexity of Algorithm 1 can be analyzed in two stages. First, the initialization stage requires the computation of the LS solution in (34), which can be handled very efficiently since the matrices C H C −1 C H and G H GG H −1 can be pre-computed offline. In particular, the complexity of the LS solution is given by O((M N + C R N + C R C ϕ )I), which results from M N I multiplications in the coherent integration and the two matrix multiplications with O(C R N I) and O(C R IC ϕ ) operations. Second, for the iterations stage, the coherent integration term has already been computed in (34) and does not change during the iterations. Similarly, the matrices C H C and GG H can be stored offline before executing the iterations. Hence, the periteration complexity of Algorithm 1 is dictated by the two matrix multiplications involved in To gain further insights, let the total number of range-angle grid points be denoted by C = C R C ϕ and assume C R = C ϕ = I for simplicity. Then, the overall complexity becomes where N iter is the number of iterations. In practice, we have observed that the algorithm converges in few tens of iterations for the RF measurement data. Consequently, the complexity of Algorithm 1 is linear in the total number of transmit symbols M N I and sub-quadratic in the number of grid points C.

IV. TRACKING-BASED DYNAMIC MAPPING
In this section, we address the actual mapping task based on a sequence of sensing measurements (range-angle charts and the corresponding identified targets) obtained by a moving UE at locations p UE l , l = 0, . . . , L − 1. Specifically, we propose and formulate a novel state model for mapping both diffuse scattering and specular reflections, which allows efficient tracking of individual scatterers or reflecting surfaces over time using IMM extended Kalman filter (EKF), where the ISTA outputs are considered as measurements. To remove second-order reflections or other second-order channel interactions, a novel measurement selection process is also presented. The tracking-based mapping allows to efficiently utilize the correlations in the sequences of measurements obtained by moving sensors, while also enables tracking of the possible changes in the radio environment.

A. Fundamentals and Rationale
At each sensing instant, the environment can be represented as a collection of discrete scatterers, which correspond to the physical points of interaction of the radio waves with the environment. These include specular reflections of the signal off the walls and other flat surfaces, and diffuse scattering off the multifaceted objects or objects with rough surfaces [52]. Tracking time evolution of the interaction points' positions can help to reconstruct the geometry of the considered environment [53]. To this end, we consider a sequence of known UE sensing positions denoted as p UE where at position l the K l points of interaction with absolute 2D coordinates are denoted by p SC l,k = (x SC l,k , y SC l,k ), k = 1, . . . , K l . Specular reflections and diffuse scattering lead to different evolutions of these interaction points: • Specular reflections: there is an unambiguous relationship between the incidence and reflection angles,so that the location of specular reflection point at each time l is a function of the positions of TX and RX at the UE, while the motion of the UE causes the shift of the specular interaction point. In this work, we consider the continuous white noise acceleration (CWNA) motion model [54], implying that the interaction point moves with constant velocity perturbed with a white noise process, to represent the kinematics of the moving scatterers. • Diffuse scattering: a relatively wide distribution of the reflection angles may correspond to each incident angle. Due to this dispersion of the reflected wave, the apparent location of the interaction point remains constant for several consecutive UE locations. To describe the kinematics of these stationary scatterers, we consider scatterer position to be constant, perturbed with a white noise process. In accordance with naming convention of [54, Ch. 6], we call this approach a continuous white noise velocity (CWNV) model in the continuation. However, the above classification is still incomplete as the indoor propagation environment may give rise to multi-bounce paths, which have different apparent mobility. To filter out such paths, we thus first present a measurement selection process, after which the specular reflections and scattering points will be tracked using an IMM EKF. Furthermore, the reflections from concave corners, albeit they can be formally classified as double specular reflections, are useful for mapping the environment. Their kinematics can be described by the CWNV model. We note that terminology-wise, the range and angle estimates obtained through the ISTA-based range-angle charting serve as the main inputs or measurements for the IMM EKF, while for the measurement selection purposes we also assume that the corresponding received signal strength (RSS) measurements are available. For readers' convenience, Observation function H l,k Jacobian matrix D a k,q Measurement association distance the basic measurement selection and IMM EKF notations are summarized in Table II. Specifically, we denote the range and angle estimates obtained through range-angle charting byd UE l,k andφ UE l,k , respectively.

B. Tracking Filter Measurement Selection Method
As noted above, the overall available mapping-related measurements for the k th scatterer obtained at the l th UE location p UE are the angle and distance estimates, respectively, obtained through range-angle charting as described in Section III, while RSS UE l,k is the corresponding RSS measurement. From the mapping point of view, the measurements originating from the single interactions are of a particular interest, hence we will harness the RSS to distinguish the first-order reflections from the higher-order interactions. In particular, we consider two alternative path loss model hypotheses, expressed in logarithmic scale: H 0 corresponds to single-order reflection model and H 1 corresponds to higher-order reflection model, with where β and β are power scaling parameters depending, e.g., on transmit power, carrier frequency and reflection coefficient, α is the path loss exponent (set to 2 under H 0 ), and n Hi is a design parameter. In Fig. 3, the RSS measurements RSS UE l,k are illustrated as a function of distance d UE l,k for all measurements of the ray-tracing data from the indoor scenario described in Section V. It can be seen that only a relatively small fraction of all measurements are actually single reflections and that at   Fig. 3. Illustration of the measurement selection approach building on two alternative pathloss models in (18) with ray-tracing data from the indoor scenario described in Section V. Distance threshold d th is also shown.
short distances, say d UE l,k ≤ d th , where d th is on the order of 2-3 meters, only single-order scattering occurs.
To estimate the unknown parameters [β, β , α ] and the classification of each measurement, denoted by ζ l,k ∈ {H 0 , H 1 }, we use the expectation maximization (EM) algorithm [55]. The measurement selection process is performed recursively over different measurement locations, using the pathloss parameter estimates from previous time steps as initial estimates for the following time step. We set the prior p prior (ζ l,k = H 0 ) = 0.5. As a result of the expectation step of the EM algorithm, each RSS measurement is associated with a posterior probability estimate p post (ζ l,k = H 0 ) ∈ [0, 1], which indicates the probability that the measurement is originated from the singleinteraction model. The final criterion for measurement selection of first-order reflection or scattering is given as where p sampl th and d th are design parameters, which affect the sensitivity of the measurement selection process.

C. IMM EKF Tracking and Measurement Association
In this sub-section, we present the IMM EKF based tracking of the scatterer locations. Each of the K scatterers is in general tracked separately, but we omit the scatterer index below for notational simplicity.
1) IMM models and weights: For generality, let us assume that there exist a set of (sub)models {M 1 , ..., M W }, and a set of corresponding prior probabilities µ j 0 = P {M j 0 }, j = 1, ..., W . At each time instant, the probability for each scatterer to switch from model i to model j are assumed known and denoted by p ij = P {M j l |M i l−1 }, i, j = 1, .., W . Both probability sets are IMM filter design parameters and are chosen to reflect the properties of the environment. At each time step, one obtains the initial conditions for each model using all previous state estimates. Also the probabilities of each model are updated at every time step. The IMM combined state estimate and the covariance are then calculated as a weighted mean of aposteriori EKF state estimates and the covariance matrices of all the models, calculated via the standard EKF prediction and update [56], where the weights are defined by the probabilities.
The details of the IMM processing steps and the propagation of the model probabilities can be found in [57].
In our case, W = 2, and as the corresponding IMM submodels {M 1 , M 2 } we consider two EKFs utilizing CWNV or CWNA motion models for the state dynamics. Note that prior model probabilities µ j 0 depend on the expected proportion between specular reflections and diffuse scattering and may generally differ for different environments and carrier frequencies.
The EKF state vector at time instant l, corresponding to the UE location p UE l , is generally denoted as s l . Then, depending on the considered motion model, the state vector includes either 2D position of a scatterer (the CWNV model) for tracking diffuse scattering points, or scatterer's position and velocity (the CWNA model) for tracking specular reflections, expressed respectively as The state vector for each model is propagated following the standard linear EKF state transition model with the non-linear measurement model [55], [56]. This can be expressed as  of the state-process noise is again one of the EKF design parameters. We note that the sizes of the state-transition and covariance matrices are defined by the used motion model, i.e., C M1 l,k , F M1 , Q M1 ∈ R 2×2 and C M2 l,k , F M2 , Q M2 ∈ R 4×4 . The corresponding update step is more involved due to the unknown association between scatterer positions at one time and the next. We first describe the measurement model and corresponding Jacobians needed in the EKF. Then we detail the solution to the data association problem.
2) Measurement model and Jacobian: The fundamental radar-based measurement vector m l,k = [φ UE l,k ,d UE l,k ] T for the k th scatterer obtained at the l th UE location p UE l contains the angle and distance estimates after the ISTA target detection step (see Section III). Measurement covariance R is in general a function of the accuracy of the angle and range measurements, while the measurement function is problem-specific and defined by the type of the sensing equipment and the techniques used. In this work, measurement errors stem from ISTA-based rangeangle charting and the related target selection process, while their distributions in the range and angle domains can be approximated with normal distribution. Furthermore, stemming from (4) and (5) The Jacobian matrix of the observation function, evaluated at the a-priori estimation of the scatterer locationp SC l,k , and the measurement residual r l,k = m l,k − h(p SC l,k ) can be shown to read where d TX l,k = p TX l − p SC l,k and d RX l,k = p RX l − p SC l,k .

3) Measurement association:
The association of scatterers observed at two consecutive time instants is unknown and must be inferred during processing. Due to the measurement noise, complicated system geometry and occasional missclassified double interactions, solving this association problem is challenging. We consider relatively slow motion of the radar or frequent measurement updates so that an apparent position of a reflection point is not expected to change dramatically during one time step, allowing us to use a proximity argument to solve this problem. Specifically, assume that {p l−1|l,k }, k = 1, ..., K l−1 , are the IMM predictions of the positions at time instant l for all scatterers that have been identified at time instant l − 1. Furthermore, assume that the coordinates of K l scatterers have been coarsely estimated from the observed angles and distances at a time instant l, using the method described in Section III, as To then associate the newly measured scatterers at step l with those tracked at the previous step l − 1 and to decide whether to continue or discontinue tracking, an association metric is calculated for each pair based on Mahalanobis distance as D a k,q = (p l,k −p l−1|l,q ) TĈ −1 l,k (p l,k −p l−1|l,q ) with k = 1, ..., K l , q = 1, ..., K l−1 .
These association distances are then used as input to a linear assignment problem to find the best global association. At the second step in scatterer association process, for each pair (p l,k ,p l−1|l,q ) under the best association, we check ifp l,k is within the p-percent confidence ellipse of IMM weighted position predictionp l−1|l,q with covarianceĈ l,k . We continue the tracking process only for associated scatterers that fulfill this condition and all others are dropped. For the newly identified scatterers (i.e., not matched to a previously detected scatter point), we initialize the IMM EKF filters. For missed scatterers (i.e., not matched to anyp l,k ), we use the prior position estimate for the next time instant.
For readers' convenience, an overall flow-chart illustrating and summarizing all the processing steps from the physicallayer IQ signals to the IMM EKF output is shown in Fig. 4.

D. IMM Smoothing
In order to utilize all the available radar-based measurements for building a map of the underlying environment, we consider an optional smoothing step. From a variety of available IMM smoothing algorithms we have chosen the approach proposed by [57], tailored to our problem, as it is relatively straightforward and does not require storing of the radar distance and angle measurements. In this approach, the IMM smoothing mimics the behaviour of the forward IMM propagation, with each of the interacting IMM sub-models using the so-called Rauch-Tung-Striebel smoothing recursion [58] combined with the model interaction. This becomes possible due to approximation of the backwards model transition probabilities by the forward model transition probabilities, stemming from the Markov properties of the model transition, as well as approximation of the smoothed probability density with the Gaussian mixture of W model conditioned smoothing densities. Such a mixture can be considered a sufficient statistics of the measurements.
As the very last stage, a thresholding operation is imposed on the refined covariances of the smoothed scatterer position estimates in order to remove the unreliable estimates. We note that the covariance of the scatterer's position is related to the SNR and to the accuracy of the angle and distance measurements. Hence, this approach allows us to use all available measurement information, automatically discarding the unreliable position estimates. This will be concretely illustrated in Section VI through the processing of both raytracing based and the actual RF measurement data.

V. EVALUATION SCENARIO AND ENVIRONMENTS
In this section, the evaluation scenario, ray-tracing environment as well as the experimental RF measurement environment and equipment, used for the validation of the proposed sensing and mapping algorithms are described. It is also noted that the complete I/Q measurement data is openly available at https://doi.org/10.5281/zenodo.4475160 [59].

A. Scenario Description
The evaluation scenario is an indoor office environment at the Hervanta Campus of Tampere University, Finland, as shown in Fig. 5(a). The considered environment consist of a corridor of 2 m wide and 60 m long with different office rooms on both sides as illustrated by the line art overlaid to all the following mapping results.
The environment is sensed along straight trajectories, one of which being shown in Fig. 5(a). The sensing related measurements are conducted with distance step of 0.5 m. In addition, these positions are deployed in both directions along the corridor, providing a total of two sets of measurements. In Fig. 5(a), the considered measurement locations as well as the most significant targets from the radar perspective are shown. We can highlight three walls of adjacent corridors that are perpendicular to the system trajectory -marked with A, B and C -located at the left side of the figure. Moreover, at the right side, three metal lockers -marked with D, E and Fare expected to be the main targets due to their notable RCS.

B. 5G NR Waveform
In both the RF measurements and the ray-tracing simulations, OFDM-based NR uplink waveform is used with the widest available channel bandwidth at the mm-wave frequency range according to [21]. Therefore, a channel bandwidth of 400 MHz with subcarrier spacing of 120 kHz is utilized. In particular, for each sensing location and scanning direction, we consider an uplink NR frequency-domain resource grid with N = 3168 active subcarriers and M = 28 OFDM symbols corresponding to an observation window of around 0.25 ms. In this case, the M consecutive OFDM symbols are coherently combined to improve the SNR of the obtained range-angle charts as described in (34). According to [60], the considered transmit waveform provides a basic radar range resolution of about 40 cm which is effectively improved by the angular measurements.
It is also shortly noted that some recent studies, such as [61]- [64], have raised the idea of joint waveform optimization to improve the sensing performance in JCAS type of systems. In this work, however, we deliberately use 3GPP 5G NR standard compliant uplink waveform with all physical channels and signal structures involved, to reflect the true waveform of NR UEs as accurately as possible.

C. Ray-Tracing Environment
For validation purposes, the RF measurement campaign results are corroborated by realistic ray-tracing simulations using Wireless Insite® [65] software. In these simulations, reproducing the scenario shown in Fig. 5(a), the measurement device follows a similar trajectory as in the RF measurement campaign with 31 test locations 1 m apart as shown in Fig. 7, using an angular scanning range from −180 • to 180 • . The TX and RX array operation is simulated through a directive  beam pattern with 3 dB beam-width of 17 • , similar to the real directive antenna systems used in the RF measurements. In addition, the same antenna separation and height are used, and the carrier frequency is 28 GHz.
The ray-tracing engine is configured to consider a maximum of 15 rays per simulation, and the number of allowed interactions is limited to six reflections and one diffraction. The walls, floor and ceiling are built using the frequency-specific materials, namely, ITU layered dry wall and floor or ceiling board. The Lambertian diffuse scattering model is applied with a scattering factor of 0.2 and a cross-polarization fraction of 0.4 to all building materials except for glass. This way, the diffuse scattering compensates for the reasonable simplifications in the environment modeling, compared to the true physical environment shown in Fig. 5(a), that were allowed to reduce the computational complexity and simulation time. However, like the results will also illustrate, the ray-tracing environment models the physical environment accurately.

D. RF Measurement Equipment
In the actual RF measurements, we use state-of-the-art mm-wave equipment to emulate the UE operation at the 28 GHz band, with selected equipment shown in Fig. 5(b). The baseline hardware platform in the measurement setup is the NI vector signal transceiver (PXIe-5840) which implements the To emulate the UE's phased array beam-steering operation, two directive horn antennas (PE9851A-20) are utilized for the TX and RX sides. These antennas are mounted on mechanical steering systems which enable to steer and direct the horns in the whole azimuth plane very accurately. According to the specifications, both horns provide a nominal gain of 20 dBi with a 3 dB beam width of 17 • . The antennas are placed at one meter above the floor level with a separation of 60 cm in order to avoid larger mutual coupling between TX and RX chains. In the TX side, two external power amplifiers (PAs) are also used that together with the antenna system facilitate an EIRP of around +20 dBm. In the RF measurements, 61 observation locations are considered as shown in Fig. 9 with an angular scanning range from −50 • to 50 • due to hardware limitations.

VI. RESULTS
In this section, we provide the ray-tracing and RF measurement based results. We start with a short example of rangeangle processing, utilizing the RF measurement data, while then put most of the focus on the actual mapping results, covering both the ray-tracing data and the RF measurement as a reference while deploys then additional smoothing and thresholding on top of them.

A. Example Range-Angle Processing Results
We start by demonstrating the proposed range-angle processing technique described in Section III using the 28 GHz RF measurement data. Specifically, the 10 th location shown in Fig. 9(e) is considered as a representative example. First, Fig. 6(a) illustrates the resulting range-angle chart when applying the considered LS estimation approach. In this case, distances up to 30 m and azimuth angles between −50 • to 50 • are investigated, using a total of C R = 391 and C ϕ = 221 range and angle cells, respectively. In addition, a Hamming window is used to improve the side-lobe suppression in the range domain. As it can be observed, the 400-MHz bandwidth used in the measurements facilitates high-accuracy range resolution, that enables to distinguish targets with mutual distances down to around 0.4 m. However, we also identify significant side-lobes in the angular-domain due to the TX and RX horn antenna patterns that degrade the target separation in this domain.
Then, a sparse representation of this range-angle chart is obtained by applying the proposed ISTA approach with regularization parameter of λ = 0.08, summarized in Algorithm 1, to better facilitate the subsequent mapping phases. The corresponding results are shown in Fig. 6(b), including also the actual detected targets that are used as input in the tracking-based dynamic mapping. In this case, a minimum target separation of 2.3 m and 20 • was considered for the range and angle domains, respectively, with a maximum of 10 detected targets per location and considering an additional dynamic range related threshold of −60 dB. As it can be seen, only the most significant target reflections remain in the sparse chart, suppressing the possible side-lobes in both range and angular domains. Hence, the ISTA approach is deployed as the main range-angle chart processing engine in the forthcoming mapping results.

B. Ray-Tracing based Mapping Results
The considered mapping methods are next assessed and validated with the ray-tracing based evaluations. As already noted, the ray-tracing environment resembles the true physical environment and the corresponding RF measurements as closely as possible -with an exception that we adopt a wider angular scanning range from −180 • to 180 • in order to seek for the maximum mapping performance. The obtained results with the different mapping methods are presented and illustrated in Fig. 7, including also the building floor plan for reference.
First, the ISTA-based static mapping approach is shown in Fig. 7(a) and compared with the grid-based static mapping of Fig. 7(b). As it can be observed, the ISTA-based approach is able to recover the most significant targets, providing a similar performance to the grid-based static method. Specifically, Fig. 7(b) illustrates the final grid-based map after applying the different averaging, filtering and thresholding stages described in [66] with the ISTA-based range-angle charts as the input. In the grid-based method, we pursue a 2D map of the simulated corridor using square cells with size of 0.2 × 0.2 m 2 , which is consistent with the range resolution of about 0.4 m, to create the initial average map. Next, as described in [66], a Gaussian kernel matrix with parameters U = V = 5 and σ = 1 is applied (see Eqs. (7)- (9) in [66]) to smooth the map and reduce the effects of noise. Finally, a thresholding stage is deployed to emphasize the most relevant targets of the environment. As it can be observed, both methods are able to fairly accurately sense the indoor scenario, providing a mapping reconstruction that clearly reflects the true layout.
Next, the same ray-tracing scenario and data is processed with the proposed tracking-based dynamic mapping approach, described in Section IV, while considering the same ISTA target detection parameters as in the previous subsection. The corresponding results are shown in Fig. 7(c) and (d), without and with IMM smoothing, respectively. These results consider the IMM-based method which deploys two separate EKF filters using the CWNV and CWNA motion models for stationary and moving scatterers, respectively. For the measurement selection procedure, we have chosen the single interaction model probability of p prior (ζ l,k = H 0 ) = 0.5 and threshold values of p sample th = 0.01 and d th = 2.3 m. In addition, the standard deviations for the single-reflection model and the multi-reflection model are defined as σ H0 = 2 and σ H1 = 10. For the initial state covariance in the EKF models, we use a variance of 10 4 m 2 for the target x and y coordinates, as well as 0.5 m 2 /s 2 for the target velocities in x and y directions when considering the CWNA model. Moreover, we choose to emphasize tracking of single-order specular reflections, and thus we initialize the target velocities at the second time step by projecting the known UE movement with a measured target angle to the movement of an assumed single-order specular reflection point at a wall. The used power spectral densities of the state covariance matrices Q c of the CWNV and CWNA models are 10 −6 m 2 /s and 0.05 m 2 /s 3 , respectively. The range and angle related measurement standard deviations are set to 10 • and 0.2 m, respectively. The initial probabilities of the scatterer motion models (µ j 0 ) for the CWNA and CWNV are set as 0.15 and 0.85, respectively.
In the forward IMM filter pass, the scatterer position estimates are tracked together with their covariance matrices. At each step, they are updated using the newly acquired angle and distance measurements. To this end, Fig. 7(c)   positions corresponding to the same scatterer cross-identified between multiple time instants are interconnected forming trajectories that are shown to follow the walls of the corridor with notable accuracy. Note that the covariance matrices of the scatterer positions become smaller as more measurement are available, hence, the position estimates in the beginning of tracking are not as reliable as those closer to the end. Also, we can still observe a few apparent scatterers beyond the walls, corresponding, e.g., to double reflections that have passed the measurement selection stage.
After the UE has finished the tracking process, it provides the tracked targets for post-processing, i.e., to the IMM smoother that is deployed to further improve the final map representation. During the backwards IMM smoothing pass, the scatterers' position and covariance estimates are updated taking into account all measurements available at the last tracking time instant. In the final map processing step, we use the smoothed covariance of each position estimate as a measure of reliability, allowing thus to extract only the most reliable ones. The corresponding final map shown in Fig. 7(d), demonstrates that a number of scatterers did not pass to the final stage due to their large covariance, especially those in the beginning of the tracking process, those far away from the TX and those corresponding to the multiple reflections.
Utilization of ray-tracing data enables extraction of groundtruth information, including the interaction coordinates of single-bounce radio paths at different walls. This allows quantitative performance evaluation of the studied methods by comparing the estimated scatterer positions with the available ground-truth data. In Fig. 8, the performance of ISTA-based measurements (including the ISTA target detection), forward IMM tracking, and backwards IMM smoothing are shown in terms of GOSPA metric [67], which considers the achieved estimation accuracy as well as the number of missed detections and false detections. The proposed IMM smoother delivers the best performance by providing almost 50% lower GOSPA compared to the ISTA-based measurements. The observed fluctuation in GOSPA over different time steps originates from the variation of measurable walls, as the UE moves along the corridor. Whenever measurements from a new wall become available, the IMM filters need to be initialized. Especially with short walls, this is a challenge as the filter states might not have enough time to properly converge. Nonetheless, the results show that despite these potential challenges, the proposed IMM approach is able to considerably improve the mapping performance compared to the original ISTA measurements.

C. RF Measurement based Mapping Results
Finally, the proposed methods are assessed and tested with RF measurements to validate the considered sensing and mapping functionality with real-world equipment and physical environment. We address how the mapping system performance is subject to the UE orientation by showing mapping results for two main UE trajectories separately. To this end, Fig. 9(a), (c) and (e) compare the grid-based static and the tracking-based dynamic mapping results when the measurement equipment is moving along the corridor from left to right, while Fig. 9(b), (d) and (f) show the corresponding mapping results for the opposite moving direction, from right to left.
It can be seen how the proposed dynamic mapping method accurately reconstructs the complex physical environment despite the fairly limited scanning range of −50 • to 50 • of the available antenna systems, compared to the ray-tracing scenario, and despite the numerous real-world effects and impairments that are present in the measurement data. We especially highlight the locations of the main corridor walls (marked with A, B and C) and the metallic lockers (marked with D, E and F) locations, both in Fig. 5 and in Fig. 9, for better interpretation of the results. We can observe how the tracking-based dynamic approach is able to accurately track the main targets of the environment, and subsequently, improve the map quality -especially when the IMM smoothing stage is deployed.
Finally, we showcase that by increasing the scanning range of the measurement setup, a more rich and further accurate representation of the environment can be obtained. In this regard, Figure 10 illustrates the smoothed tracking-based dynamic mapping results when a wider scanning range of −90 • to 90 • is used in the measurements. As it can be observed, the wider scanning range provides a more complete representation of the indoor map in comparison with the map shown in Fig. 9(f). Overall, the mapping results with real-world RF measurement data demonstrate the applicability of the proposed methods also in true complex physical environments, allowing to extract situational awareness in an efficient manner. In general, when the future mobile networks further evolve towards the 100 GHz bands, and perhaps even beyond, the relative role of diffuse scattering is likely to increase [24]. The novel state formulation described in this article -taking into account both the specular and diffuse componentscombined with the IMM filtering framework is offering a versatile tool for high-accuracy radio environment mapping also in such networks. Demonstrating that through concrete RF measurements is one important ingredient in our future work. VII. CONCLUSIONS This paper investigated the sensing and environment mapping prospects of mm-wave 5G NR and beyond networks with specific emphasis on the UE side for mobile mapping applications. In the considered framework, the UE operates as a joint communication and sensing system, scanning its surrounding environment while steering its beam pattern towards different directions and observing the target reflections. First, we derived a novel LS-based processing technique to estimate sparse rangeangle charts to facilitate the subsequent mapping processing. Then, an advanced dynamic mapping approach was presented, building on novel state model with diffuse and specular scattering together with IMM-EKF filtering and smoothing. While the sparse range-angle charts already allow for baseline static mapping capability, the dynamic tracking-based approach that builds on the IMM-EKF processing solution incorporating both specular reflections and diffuse scattering allows to track individual scatterers over time. This together with the IMM smoothing and the described novel scatterer measurement selection and association methods provide an efficient framework to reconstruct and map complex physical environments. The applicability of the proposed methods and algorithms was then assessed and evaluated, through both ray-tracing simulations and actual RF measurements at the 28 GHz band in practical indoor office type of an environment. The obtained results demonstrate the good applicability of the proposed methods, while the IMM-EKF based dynamic tracking solution was shown to clearly outperform the more simple static approaches. Our future work will consider extensions to 3D sensing and mapping, while extending also to the actual SLAM where the coordinates of the sensing device are also unknown. Additionally, the future work contains generalizing the presented methods for multiple simultaneous transmit and/or receive beams, and carrying out measurements at 60-100 GHz bands.

APPENDIX A LS SOLUTION
The LS approach to estimate the range-angle chart B in (12) is formulated as After vectorization, we have where x m vec (X m ), y m vec (Y m ), b vec (B) and vec (·) denotes the vectorization. Further simplifying (27) yields where is a known matrix. Defining then y = y T 0 , . . . , y T M −1 The LS estimate of the range-angle chart in (32) can thus be obtained as where (·) † denotes matrix pseudo-inverse. In the regime of large number of OFDM subcarriers/symbols, the LS solution can be approximated as (see Appendix B for the derivations): which consists of coherent integration over the M received symbols, range compression/matched filtering and antenna pattern correlation operations.

APPENDIX B APPROXIMATED LS SOLUTION
In this Appendix, we show how (33) can be approximated as (34) using the Kronecker structure in (29). Suppose Φ has full column rank -an assumption that is well-justified due to randomness of data symbols x m and large number of OFDM subcarriers/symbols so that N IM > C R C ϕ . Then, we can rewrite (33) as Going from (35) to (36), we use the law of large numbers for sufficiently large M by assuming that the constellation has an average magnitude of 1. The expression in (36) can be further re-written as which finally yields As seen from (39), the LS solution corresponds to a matched filtering operation. In particular, the observation matrix is first multiplied by the conjugate of the OFDM transmit symbols to cancel out their effect. Then, the resulting matrices are coherently integrated over M symbols relying on the assumption of negligible Doppler over the duration of M symbols. Next, left multiplication by C H C −1 C H represents range compression/matched filtering via the rangedependent frequency-domain steering matrix C. Notice that the columns of C coincide with those of an N × N DFT matrix. This means that matched filtering via C H C −1 C H is essentially a normalized IDFT operation over the columns of X * m Y m , which is a standard approach in OFDM radars [42], [43]. Similarly, right multiplication by G H GG H −1 in (39) represents correlation/matched filtering with the antenna pattern. Hence, (39) can be interpreted as matched filtering in the range-angle domain, which essentially maximizes the SNR of the target detection.