METHOD AND SYSTEM FOR THE SPARSE RECONSTRUCTION OF THE MICRO-DOPPLER SPECTRUM IN JOINT COMMUNICATION AND SENSING APPLICATIONS

Abstract

The present invention refers to a method and system for joint communication and reconstruction of the micro-doppler time-frequency spectrum from sparse channel measurements.

Claims

1. A computer implemented method for joint communication and reconstruction of a micro-Doppler time-frequency spectrum from sparse channel measurements, wherein wireless communication signals, including channel estimation fields, are transmitted through a multi-path channel, and the reflections or refractions of the transmitted signal are received, comprising: i) estimating the channel impulse response (CIR) of the multi-path channel, wherein the CIR contains complex channel gains for each path of the multi-path channel, to obtain a plurality of CIR estimates corresponding to irregularly spaced CIR values; ii) resampling the available channel impulse response (CIR) estimates to obtain an incomplete regular grid using a resampling technique, the incomplete regular grid comprising CIR samples regularly spaced in time with possibly some missing samples; and iii) performing a sparse reconstruction of a Fourier transform of the incomplete regular grid in the time domain for reconstructing the micro-Doppler time-frequency spectrum.

2. The method according to claim 1, wherein said resampling is performed on a slotted sliding window with slot duration (T) and window length slots (W), and wherein said steps i), ii) and iii) are repeated for each subsequent window.

3. The method according to claim 2, wherein said time slot duration (T) is selected as T=c/(4f.sub.ov.sub.max) where v.sub.max is the desired maximum micro-Doppler velocity resolution in the spectrum, c is the light speed and ft is the carrier frequency.

4. The method according to claim 2, wherein the following steps are performed between said estimating and resampling steps: a) setting a threshold number of CIR measurements-per-window, needed to reach the desired micro-Doppler reconstruction quality; and b) transmitting a number of additional CIR estimation fields to meet the threshold number of measurements.

5. The method according to claim 4, further comprising a step of scheduling the additional CIR estimation fields to be transmitted in the current window according to a predefined scheduling policy.

6. The method according to claim 5, wherein said scheduling policy is to transmit K additional CIR estimation fields in the last K slots of the time window.

7. The method according to claim 5, wherein the scheduling is performed at half of the window duration (W/2) and the subsequent window is shifted forward by W/2 slots.

8. The method according to claim 5, wherein, after the scheduling has been performed, for any CIR estimate extracted from a communication packet that is received after the scheduling operation, the first scheduled CIR estimation field is removed from the schedule.

9. The method according to claim 2, wherein said resampling step comprises selecting, for each slot, the CIR value sampled at the time instant closest to the slot center and, if no sample is obtained in a slot, the CIR window slot sample is considered missing.

10. The method according to claim 1, wherein said step of performing sparse reconstruction is performed separately for each signal propagation path.

11. The method according to claim 2, implementing a reconstruction algorithm comprising: a) building an W by W inverse Fourier basis matrix (B) wherein W is the number of slots in the time window; b) building a reduced inverse Fourier matrix (F) containing the rows of B with indices corresponding to the non-missing CIR samples in the window; c) posing an optimization problem such that its solution, a vector H of dimension W, is the Fourier transform of the complete CIR measurement window, and such that it enforces the sparsity of H; and d) solving said optimization problem to obtain H and computing the spectrum as H.sup.2.

12. The method according to claim 11, wherein the algorithm used for solving the optimization problem uses the Iterative Hard Thresholding (IHT) method.

13. The method according to claim 1, wherein said reconstruction is performed only on the path yielding the highest received power.

14. The method according to claim 1, wherein said reconstruction is performed on a subset of the paths contained in the estimated CIR, being the subset containing the contribution of a target of interest.

15. The method according to claim 14, wherein the spectra from the different paths of the subset are combined by summing the micro-Doppler spectra obtained from the Fourier transforms.

16. A system for joint communication and reconstruction of the micro-Doppler time-frequency spectrum from sparse channel measurements, comprising: a) a transmitter configured for transmitting wireless communication signals through a multipath channel, including channel estimation fields; b) a receiver, configured for receiving a wireless signal which is the reflection or refraction of the transmitted signal; and c) a processor, configured for implementing a method according to claim 1.

17. The system according to claim 16, wherein said transmitter and said receiver are a single transceiver sharing a same antenna array working in full-duplex mode.

18. The system according to claim 16, wherein the estimated CIR is from a backscatter channel.

19. The system according to claim 16, wherein the transmitter is equipped with an antenna array and phase shifters for directional beamforming.

20. The system according to claim 19, wherein the CIR is estimated for each of the different beampatterns used during the transmission.

21. (canceled)

22. (canceled)

Description

BRIEF DESCRIPTION OF THE FIGURES

[0029] Hereinafter in this description, reference will be made to the drawings shown in the enclosed figures, wherein:

[0030] FIG. 1 is a comparison between the traditional CIR-based human sensing and the present invention.

[0031] FIG. 2 is a visual representation of the D spectrum computed using the present invention on two different CIR paths, one containing the reflection from a person's torso, the other capturing the D signature of the leg. The total D is obtained summing together these contributions.

[0032] FIG. 3 shows an example injection procedure with M.sub.s=8, W=16. Four sensing units are scheduled after P1 and P2. Then, as three reflected communication packets are received, they are reused and the first three scheduled sensing units are not transmitted. The fourth sensing unit is instead injected in the last slot.

[0033] FIG. 4 shows an exemplary block diagram of a system according to the invention.

[0034] FIG. 5 shows walking D spectrograms and RMSE for different levels of sparsity, obtained by uniformly removing samples for each window.

[0035] FIG. 6 shows exemplary traffic patterns from the pdx/vwave dataset.

[0036] FIG. 7 shows per-class F1 scores obtained for different values of M.sub.s and RAPID on the tested dataset.

[0037] FIG. 8A shows overhead of the invention for different values of M, in the three traces of the pdx/vwave dataset.

[0038] FIG. 8B shows overhead vs. average HAR F1 score for different values of M.sub.s.

[0039] FIG. 9A shows per-class F1 scores aggregating a different number of paths Q.

[0040] FIG. 9B shows per-class F1 scores changing , the IHT sparsity parameter.

TECHNICAL BACKGROUND

[0041] In the following, a brief description of the CIR model for communication systems used for sensing is provided. Then a baseline approach is described that allows to track the movement of people in the environment and extract their D signatures using regularly sampled CIR information. This forms the basis of the present invention, which eliminates the requirement of fixed Inter-Frame Spacing (IFS) and enables ultra low-overhead ISAC.

Sensing in Communication Systems

[0042] Capturing the movement features of humans in the environment requires an analysis of the reflections of the transmitted signal from their bodies, which is usually carried out applying signal processing techniques to the CIR. In the case of standard communication systems, for example, network devices implementing the IEEE 802.11 standard at 2.4 or 5 GHz, communication is typically performed omni-directionally. For this reason, and due to the intrinsic lower distance resolution of such frequencies, super-resolution techniques would be needed to accurately localize and sense the subjects. In the case of mmWave systems instead, due to the high path loss occurring at mmWave frequencies, directional communication is employed by means of transmitter and receiver beamforming, typically by means of phased antenna arrays. The transmitter and the receiver use suitable BP configurations of their antenna arrays to maximize the signal strength. To successfully sense with a mmWave system, at least one of the BPs has to illuminate the subjects, as only in this case the reflected signal carries detectable information about the movement signature. To this end, in the testing of the present invention a setup is considered where an AP transmits packets and is able to collect the reflections of its own signal, after being reflected by objects (including humans). This reflection is preferably collected by the receive array of the AP itself using a quasi-omnidirectional BP. This requires full-duplex capabilities, as is common in ISAC scenarios, which in the simplest form can be achieved with a Multiple Input Multiple Output (MIMO) system in a mixed configuration with one RF chain as transmitter and another as receiver. The full-duplex configuration is not critical for the functioning of the proposed method and only constitutes a preferred embodiment used in the evaluation of the invention. The sensing-oriented CIR estimation fields, which are denoted by sensing units, can either be transmitted independently (injected) or piggybacked by appending them as a trailer to the Physical Layer (PHY) communication packets. mmWave standards implement beam training mechanisms that help to establish a communication link by testing different BP combinations and then selecting the best one. Such functionality is supported by all mmWave standards. For example, 5G-NR, use Synchronization Signal Block (SSB) and CSI-RS for beam management, while WLAN systems adopting the IEEE 802.11ad/ay standards use channel estimation and training fields (CEF and TRN, respectively) to obtain accurate CIR information. The framework to extract sensing information from CIR measurements can be applied regardless of the specifics of the standards.

CIR Model

[0043] Channel measurements in communication systems contain information about the environment. Depending on the communication system, sensing could be performed using, e.g., the 5G-NR and IEEE 802.11n/ac/ax Orthogonal Frequency Division Multiplexing (OFDM) Channel State Information (CSI), which contains the channel gains for each OFDM subcarrier, or the IEEE 802.11 ad/ay Single Carrier (SC) CIR.

[0044] Both communication schemes are suitable for human sensing: (i) in OFDM systems, the base stations can send frequent downlink packets to estimate the channel, possibly using different BPs if transmit beamforming is enabled, while (ii) in SC systems, such as IEEE 802.11ad/ay, in-packet beam tracking can be used, so that specific fields called training fields (TRN), each using a different BP, can be appended to communication packets. In the following, the focus will be on SC CIR, and it will be shown how to extract the D effect of human movement. However, previous works have demonstrated that similar processing can be performed with OFDM CSI, and the present invention is general enough to be applied in both cases. In addition, the CIR model presented in the following assumes that the transmitter device can implement transmitter beamforming to direct the signal energy along preferred directions. This is coherent with the preferred implementation of the system but the proposed method is general enough to be applied to systems that do not implement transmitter beamforming.

[0045] The CIR contains the complex channel gains for each path l=0, . . . , L1 along which the signal travels. Each path is associated to a specific distance from the AP, according to the relation d.sub.l=c.sub.l/2, with .sub.l being the delay of the reflection from path fand c the speed of light. The ranging resolution of the system, i.e., its capability to resolve the distance of the reflectors causing different signal paths, is given by d=c/2B, termed range bin, where B is the bandwidth of the transmitted signal. Moreover, the CIR depends on the specific BP used during the transmission, denoted by b=0, . . . , N.sub.BP1.

[0046] For carrier frequency f.sub.o, the CIR along , using BP b at time t is

[00001]$\begin{array}{cc}{h}_{,b}\left(t\right)={.Math.}_{p=1}^{P_{}\left(t\right)}{a}_{,b}^p\left(t\right)\exp \left\{-j2\frac{2f_o}{c}\left[d_{}+{}_0^t{v}_{}^p\left(x\right)\mathrm{dx}\right]\right\}.& \left(1\right)\end{array}$

[0047] In Eq. (1), (t) is the number of reflectors that are located within d from , whose contributions are overlapping in path , while

[00002]$v_{}^p$

is the radial velocity (by convention,

[00003]$v_{}^p$

has a positive sign when the reflector moves away from the AP) of the p-th reflector. The quantity

[00004]$a_{,b}^p\left(t\right)$

contains the complex gain due to the joint effect of the transmitter BP, the reflectivity of the object and the signal attenuation.

Micro-Doppler Extraction

[0048] The extraction of the D spectrum from multiple, concurrently moving subjects requires tracking the position of each person in the physical space, in order to separate their individual contributions to the CIR. Then, a spectral analysis over different CIR samples yields the desired D signature.

People Tracking.

[0049] In a possible embodiment of the present invention, a people tracking phase can precede the micro-Doppler extraction, in order to perform the extraction separately for each subject to be sensed. People tracking is performed by extracting measurements of each person's distance and angular position with respect to the AP across time. This process consists of (i) removing the background contribution to the CIR by subtracting the average CIR across a suitable time interval, (ii) selecting the locally strongest reflection paths in the CIR (peaks) and obtaining the corresponding distance d, (iii) computing the Angle of Arrival (AoA), , of the reflection from the correlation between the different BPs gains and the strengths of the CIRs across the whole angular Field-of-View (FoV). The approach in (iii) requires the BP shapes to be estimated in advance, and is based on the intuition that different BP illuminate different possible reflectors in the environment, depending on their position, so that one can expect to receive the reflected signal from a certain subject only when a BP pointing in his/her direction is used. The resulting set of distances and angles represent candidate positions of humans in the environment. A multi-target tracking method such as a Joint Probabilistic Data Association Filter (JPDAF) allows smoothing the trajectories of the subjects rejecting noise and clutter.

[0050] The tracking phase provides an estimate of the position of each subject, at every time instant t, which is denoted by [{circumflex over (d)}(t), {circumflex over ()}(t)] using the superscript {circumflex over ()} to differentiate between the estimate and the true position. Due to step (ii) above, which selects the local peaks in the received power, typically only the reflections from the torso can be tracked when a person moves in the environment. This is because the head and the limbs cause much weaker reflections, whose contributions can only be detected in the D spectrum. Nevertheless, tracking the spatial location of the torso is crucial for the D extraction as it allows separating the contributions of the multiple subjects in the environment. Once {circumflex over (d)}(t) and {circumflex over ()}(t) are determined, they are used to select the path * and the BP b* that correspond to the distance and angular position of the subject, respectively.

[0051] Then, the CIR waveform that contains the D effect of the person's movement is (t). Besides enabling the separation of multiple subjects, this operation makes mmWave sensing systems much more robust to changes in the environment than sub-6 GHz devices. While the latter are heavily affected by second-order reflections on walls and objects, mmWave sensing mostly relies on the line-of-sight path between the person and the transmitter/receiver, with little contribution from the external environment.

Micro-Doppler Spectrum.

[0052] Human movement causes a small scale Doppler effect on the reflected signal due to the different body parts, which possess different velocities and follow different trajectories. This is referred to as D effect and causes a measurable frequency modulation on the reflection. High frequency signals such as mmWave communications are particularly affected by the D modulation due to their small wavelengths. Various techniques from time-frequency analysis can be applied to analyze the D, obtaining spectrograms showing the time evolution of the signal energy contained in the different frequency bands of interest. The most popular and computationally efficient of such methods is the Short Time Fourier Transform (STFT) of (t), which consists in applying a windowed Fourier Transform (FT) on partially overlapping portions of the CIR. STFT processing needs to be applied on a window of W subsequent estimates of the CIR, computed with a fixed sampling period of T, seconds, provided that the time spanned by the window is short enough to consider the movement velocity of the reflectors constant for its whole duration. Note that this operation allows to detect and separate the velocities of the (t) reflectors, whose contributions are overlapping in path when considering a single estimate of the CIR. The choice of T.sub.c impacts the frequency resolution of the STFT, f.sup.d=1/(WT.sub.c), and its maximum measurable frequency,

[00005]$f_{m\mathrm{ax}}^d=1/\left(2T_c\right).$

Using the relationship between the Doppler frequency and the corresponding velocity, one can obtain the velocity resolution and the maximum observable velocity as v=c/(2f.sub.oWT.sub.c) and v.sub.max=c/(4f.sub.oT.sub.c). To fully capture the range of velocities of interest for human movement, the typical approach is to select T.sub.c such that v.sub.max is sufficiently high that is covers the velocities that can occur in the human activities of interest, which may vary depending on the application.

[0053] Previous work assumes that the constraint of a fixed T.sub.e is met, which does not hold in realistic communication scenarios, where the packet transmissions are scheduled according to the needs of the communication protocols rather than sensing accuracy. Traffic patterns are typically bursty and irregular and thus cannot be used by existing methods for human sensing. Instead, dedicated time slots need to be reserved for the transmission of sensing units, which is incompatible with the random access CSMA/CA MAC commonly used in IEEE 802.11. Conversely, the present invention is the first approach that does not require any specific pattern in the transmission of the sensing units, thus enabling a true ISAC where communication packets are exploited for sensing whenever possible, and only a minimal additional overhead is introduced when necessary.

DETAILED DESCRIPTION OF PREFERRED EMBODIMENTS

[0054] The present invention will be described hereinafter by making reference to the above-mentioned figures.

[0055] The algorithm to recover the D spectrum from irregular and sparse CIR sampling patterns according to the present invention is now presented. The processing steps of the method according to the present invention, compared to traditional CIR-based sensing methods are shown in FIG. 1.

[0056] In general terms, a method according to the present invention allows joint communication and reconstruction of the micro-Doppler time-frequency spectrum from sparse channel measurements, wherein wireless communication signals, in the form of packets including channel estimation fields are transmitted through a multi-path channel, and the reflection or refraction of the transmitted signal are received.

[0057] Such method essentially comprises the following steps: [0058] i) estimating the channel impulse response (CIR) of the multi-path channel; [0059] ii) resampling the available channel impulse response (CIR) estimates to an incomplete regular grid of CIR measurements using a resampling technique; [0060] iii) performing sparse reconstruction of the Fourier transform of the incomplete regular grid in the time domain.

[0061] In particular, after CIR estimation and, eventually, people tracking, for which, for example, the standard JPDAF technique can be adopted, a resampling technique to approximate the irregularly spaced CIR values with a regular grid whose sampling interval is chosen according to the desired D resolution. Due to irregularity of the original sampling process, the approximated regular sequence may be incomplete, containing missing values that need to be filled in the subsequent processing steps.

[0062] The recovery of the D spectrum from the resampled incomplete regular grid of CIR measurements is formulated as a sparse recovery problem. For this, the intrinsic sparsity of the propagation channel is leveraged, which leads to a small number of signal reflections from the human body, each carrying information about a different body part.

[0063] According to preferred embodiments, the sparse recovery problem is preferably solved with an algorithm using the Iterative Hard Thresholding (IHT) method for each signal propagation path (CIR path), and aggregate the results to obtain the final D spectrum. The aggregation may be, for example performed by summing the micro-Doppler spectra obtained from the Fourier transforms.

[0064] According to possible embodiments of the invention, when communication traffic is absent or too scarce to obtain an accurate reconstruction, according to the invention short sensing units can be injected into the (idle) channel to overcome the problem. Thanks to the sparse reconstruction of point (2), the amount of units that need to be injected is minimal and can be tuned to trade off between overhead and sensing accuracy.

[0065] In particular the resampling is performed on a slotted sliding window with predetermined slot duration (T) and window length slots (W), and wherein the steps i), ii) and iii) are repeated for each subsequent window.

[0066] Then, the following steps are performed between said estimating and resampling steps: [0067] a) setting a threshold number of CIR measurements-per-window, needed to reach the desired micro-Doppler reconstruction quality; [0068] b) transmitting a number of additional CIR estimation fields to meet the threshold number of measurements, if the communication traffic is not sufficient.

[0069] According possible embodiments of the invention, the time slot duration T is selected as T=c/4 f.sub.ov.sub.max where v.sub.max is the desired maximum micro-Doppler velocity resolution in the spectrum, c is the light speed and f.sub.o is the carrier frequency.

[0070] According possible embodiments of the invention, the resampling step comprises selecting, for each slot, the CIR value sampled at the time instant closest to the slot center and, if no sample is obtained in a slot, the CIR window slot sample is considered missing.

[0071] According to possible embodiments, the method further comprises a step of scheduling the additional CIR estimation fields to be transmitted in the current window according to a predefined scheduling policy. For example, such scheduling policy is to transmit K additional CIR estimation fields in the last K slots of the time window. A possible implementation of the scheduling policy provides that it is performed at half of the window duration W/2 and the subsequent window is shifted forward by W/2 slots.

[0072] According to possible embodiments, once the scheduling has been performed, for any CIR estimate extracted from a communication packet that is received after the scheduling operation, the first scheduled CIR estimation field is removed from the schedule.

[0073] According to possible embodiments, the method of the invention implements a reconstruction algorithm comprising: [0074] a) building an W by W inverse Fourier basis matrix (B) wherein W is the number of slots in the time window; [0075] b) building a reduced inverse Fourier matrix (F) containing the rows of B with indices corresponding to the non-missing CIR samples in the window; [0076] c) posing an optimization problem such that its solution, a vector H of dimension W, is the Fourier transform of the complete CIR measurement window, and such that it enforces the sparsity of H; [0077] d) solving said optimization problem to obtain H and computing the spectrum as H.sup.2.
Preferably, the reconstruction can be performed only on the path yielding the highest received power or, in alternative, on a subset of the paths contained in the estimated CIR, being the subset containing the contribution of a target of interest.

[0078] In the following, a more detailed description of the features described above will be given with the support of related mathematics concept and principles the following description refers to possible combinations of the steps step of the method as above described, as a single embodiment. This is done in order to render the description as clear and fluent as possible. However, it is to be intended that different embodiments could be provided, using different combinations of features.

CIR Resampling

[0079] According to the invention the CIR are sampled at time instants t.sub.i, which coincide with the reception of the reflections from the i-th transmitted packet. To reconstruct the D spectrum from CIR samples which are randomly distributed in the time domain, it is firstly resampled the CIR to obtain regularly spaced samples with a fixed granularity T.sub.c where possible.

[0080] To do so, it is resorted to the slotted resampling technique, which allows approximating a sequence of randomly spaced samples into a regular grid with missing values. N.sub.s consecutive samples are considered, obtained at the time instants t.sub.0, t.sub.1, . . . , t.sub.N.sub.s and denote by 0, T.sub.c, 2T.sub.c, . . . , (K1)T.sub.c the regular grid with step size Tc. Slotted resampling constructs a new CIR sample sequence (kT.sub.c) where the CIR values are obtained from the original sequence (t.sub.i) as follows. Time bins (or intervals) of length T, are centered on each time instant of the regular grid, i.e., bin k is .sub.k=[kT.sub.cT.sub.c/2,kT.sub.c+T.sub.c/2), with center kT.sub.c. Then, the value of the CIR corresponding to the grid value k is either (i) selected among the values of the original sequence whose sampling times fall inside bin k, taking the one whose sampling time is the closest to the bin center, or (ii) considered as a missing value if no samples of the original sequence fall inside bin k. Specifically,

[00006]$\begin{array}{cc}{h}_{,b}\left({\mathrm{kT}}_c\right)=\{\begin{array}{cc}0& \mathrm{if}\left\{t_i.Math.t_i{}_k\right\}=,\\ h_{,b}\left(t_k\right)& \mathrm{otherwise},\end{array}& \left(2\right)\end{array}$

where the 0 values represent missing samples and

[00007]$\begin{array}{cc}{t}_k=\arg \underset{\left\{t_i.Math.t_i{}_k\right\}}{\min }.Math.{\mathrm{kT}}_c-.Math..& \left(3\right)\end{array}$

[0081] The resulting, regularly spaced sequence of CIR samples is used to reconstruct the D spectrum of the subject. However, due to the missing samples which are set to 0, a plain application of the STFT (as described above) would lead to a corrupted spectrum.

[0082] In the following the solution to this problem is detailed, which is based on sparse recovery techniques.

Sparse AD Recovery Problem Formulation

[0083] Several methods exist to tackle the problem of computing the power spectrum of non-uniformly sampled signals. Such approach belongs to the category of sparsity-based approaches, in which the sparsity of the signal in the frequency domain is leveraged to drastically reduce the number of measurements needed for an accurate reconstruction of the spectrum. Windows of length W samples (window size) every samples (overlap) from the sequence h.sub.l,b(kT.sub.c) are selected, choosing =W/2. Due to the slotted resampling process, each window may contain missing samples. The set of indices of the available samples contained in the m-th window is denoted by .sub.m.

[0084] Then, it is defined vector .sub.b(m).sup.m|, containing the available CIR samples in the m-th window, and vector .sub.b(m).sup.W, representing the complete m-th CIR window, which is only partially known due to the missing samples. It is also denoted by F.sub.inv the inverse Fourier matrix, whose element in position (g, l) is given by (F.sub.inv).sub.gl=(1/{square root over (W)})exp(j2gl/W), g,l=0, . . . , W1 while U.sub.m=[u.sub.i.sup.T], .sub.m, iUm is the matrix that selects the rows of F.sub.inv whose indices are in m. .sub.m. u.sub.i is the vector of all zeros but the i-th component, which equals 1.

[0085] The following relation holds between the incomplete CIR window, .sub.b(m), and the Fourier Transform (FT) of the full CIR window, .sub.b(m).sup.W, which it is aimed to recover in order to compute the D spectrum,

[00008]$\begin{array}{cc}{h}_{,b}\left(m\right)=U_m{\tilde{h}}_{,b}\left(m\right)=U_m{F}_{\mathrm{inv}}{H}_{,b}\left(m\right)={}_m{H}_{,b}\left(m\right),& \left(4\right)\end{array}$

wherein the last step matrix .sub.m=U.sub.mF.sub.inv is used as a shorthand notation. Given Eq. (4), the aim is to recover H.sub.l,b(m) from the incomplete measurement vector .sub.b(m), which is a typical sparse recovery or compressed sensing problem. In this framework, it has been proven that recovering the FT of the desired signal is possible if the latter is sparse in the frequency domain, i.e., the FT only contains a low fraction of non-zero elements. To verify that this sparsity assumption holds in this case, Eq. (1) can be rewritten after the resampling and windowing operations, so that the i-th sample of the complete m-th window is given by

[00009]$\begin{array}{cc}{\left[{\tilde{h}}_{,b}\left(m\right)\right]}_i={.Math.}_{p=1}^{P_{}\left(m\right)}{a}_{,b}^p\left(m\right)\exp \left\{-j4\frac{f_o}{c}\left[d_{}^p+\left(m+i\right)T_c{v}_{,m}^p\right]\right\},& \left(5\right)\end{array}$

where

[00010]$v_{,m}^p$

is the radial velocity of the p-th reflector in path during window m, and its distance from the AP. Here, it is assumed that the velocity of each reflector can be considered constant during a window. In addition, it is also considered that the reflective coefficients and the number of reflectors are constant. This is reasonable for the considered setup, where the reflectors are parts of the human body, which typically move slowly compared to the duration of a window WT.sub.c.

Single-Path Sparse Recovery

[0086] Given the model from Eq. (4), the reconstruction of the CIR FT along each path can be posed as a sparse recovery problem. Specifically, it is sought a vector .sub.b*(m) which is a solution to Eq. (4) while being as sparse as possible, coherently with the above discussion.

[0087] Considering the BP b* pointing in the direction of the target, the desired FT of .sub.b*(m) can be found as a solution of the optimization problem

[00011]$\begin{array}{cc}{H}_{,b^{*}}\left(m\right)=\arg \underset{H}{\min }{.Math.H.Math.}_{0}s.t.{.Math.h_{,b^{*}}\left(m\right)-{}_{m}H.Math.}_2,& \left(6\right)\end{array}$

where .sub.0 denotes the .sub.0-norm of a vector, i.e., the number of its non-zero components. The constant >0 can be estimated from the noise in the CIR, using a training dataset.

[0088] An approximate local solution to Eq. (6) can be found using fast greedy algorithms. It is adopted the IHT, which solves

[00012]$\begin{array}{cc}{H}_{,b^{*}}\left(m\right)=\arg \underset{H}{\min }{.Math.h_{,b^{*}}\left(m\right)-{}_{m}H.Math.}_2^{2}s.t.{.Math.H.Math.}_0,& \left(7\right)\end{array}$

where is a pre-defined sparsity level parameter. The algorithm involves an iterative gradient descent step on the quadratic term in Eq. (7), followed by a thresholding operation:

[00013]$\begin{array}{cc}{\hat{H}}^{\left(n+1\right)}{}_{}\left[{\hat{H}}^{\left(n\right)}+{}_m^T\left(h_{,b^{*}}\left(m\right)-{}_m{\hat{H}}^{\left(n\right)}\right)\right],& \left(8\right)\end{array}$

where n is the iteration index and TO is the hard-thresholding operator, which sets to 0 all the components of the argument vector except the largest ones in terms of the Euclidean norm. is a learning rate parameter which can be tuned to improve the convergence properties. The iterative process is stopped whenever .sup.(n+1).sup.(n).sub.2< or when a maximum number of iterations, nmax, is reached. According to the present invention, is a key parameter, which is strictly related to the number of reflectors P.sub.l(m): indeed, as IHT reconstructs a vector which has at most non-zero elements, is an upper bound for P.sub.l(m), and it can be thought of as the maximum number of reflectors per path that are allowed to be reconstructed. (can be tuned in order to obtain better D reconstruction. The sparse recovery algorithm is summarized in Algorithm 1.

TABLE-US-00001 Algorithm 1 - Single path sparse recovery. Input: * (m), , n.sub.max, , . Output: * (m). 1: Collect the set of available samples indices . 2: Build matrices U.sub.m = [u.sub.i.sup.T], i and 3: Compute .sub.m = U.sub.m . 4: Set .sup.(0) = 0, n = 0, .sup.(0) to any value > . 5: while n < n max or (n) > do 6: .sup.(n+1) Eq. (8) 7: .sup.(n+1) || .sup.(n+1) .sup.(n) ||.sub.2 8: n n + 1 9: end while 10: return .sup.(n)

[0089] According to the compressive sensing theory, the reconstruction performance of IHT (and in general of any recovery algorithm) degrades as the number of available measurements, |.sub.m|, decreases.

[0090] Theoretical results show that the minimum number of measurements needed to reconstruct .sub.b*(m) is ( log(W/)), although the exact number needs to be estimated empirically as it also depends on the level of noise present in the signal. In the following it will be demonstrated that the present invention can achieve excellent D reconstruction with as low as W/8 measurements per window, thanks to the high sparsity of the mmWave CIR.

Multi-Path Aggregation

[0091] The moving body of a person causes several reflections that affect more than one CIR path, as already discussed. Using the described procedure, the present invention is able to retrieve the contribution of each path .sub.b*(m) to the D. Since the different body parts contribute to the D in different paths, to fully capture human movement it is needed to combine the information from the different paths. It is called Q the number of paths considered in obtaining the D spectrum. The spectra obtained by the path caused by the torso, l*, are aggregated with the Q/2 paths preceding l* and the Q/2 subsequent paths, as they may contain the contributions of the other body parts. The expression of the total D spectrum is

[00014]$\begin{array}{cc}D\left(m\right)={.Math.}_{={}^{*}-.Math.Q/2.Math.}^{{}^{*}+.Math.Q/2.Math.}{.Math.H_{,b^{*}}\left(m\right).Math.}^2,& \left(9\right)\end{array}$

where the squared magnitude is applied element-wise. In addition normalization of the spectra in the range [0, 1] is applied by computing

[00015]$D\left(m\right)\frac{D\left(m\right)-{\min }_i{D}_i\left(m\right)}{{\max }_i{D}_i\left(m\right)-{\min }_i{D}_i\left(m\right)}.$

[0092] Note that Eq. (9) entails solving Q optimization problems of the form in Eq. (7), however, the Q problems can be parallelized as they are completely independent. Decomposing the full D spectrum reconstruction problem into Q subproblems effectively allows applying sparse recovery techniques, which in turn lead to a significant reduction of the number of measurements needed.

[0093] The value of Q is selected here due to physical considerations and validated in practice. The D vectors from Eq. (9) can be collected in sequences, one every 8 slots, forming D spectrograms of arbitrary length, depending on the specific application that is being performed, e.g., activity recognition, fall detection, gait segmentation, etc. In the following, reference will be made to the number of D vectors considered in such spectrograms as A.

Sensing Unit Injection

[0094] The present invention can exploit the sensing units in sparsely distributed communication packets to recover the D spectrum of human movement.

[0095] However, during communication between the AP and one or more terminals it may happen that the AP remains silent for longer than the duration of a processing window, W, or that the number of received packets is less than the minimum number of measurements required for an accurate D reconstruction. In these cases, the sparse recovery algorithm cannot recover H.sub.l,b(m) as the available sensing units are insufficient. To tackle this problem, the system is allowed to inject sensing units into the channel whenever the number of communication packets is not sufficient for Algorithm 1 to work, i.e., when less than two CIR measurements are available in a window. Different from existing ISAC frameworks, however, the used sparse recovery approach allows to introduce a minimal amount of overhead, as the D spectrum can be recovered from a number of CIR samples which is much lower than the full length of the window W. Note that for the injection of a sensing unit it is sufficient to transmit the necessary CIR estimation fields, without any preamble and header as used in conventional packets, since the unit is only received at the AP itself and contains a known waveform.

Basis of the Injection Algorithm

[0096] In the following, the proposed injection procedure will be presented as if both communication packets and sensing units were transmitted at times that lie on a uniform grid with spacing T.sub.c. This simplification is based on the fact that the already described slotted resampling process is used. Therefore, the injection process can be described in terms of windows of size W, where each value in the window occupies a slot which is a multiple of T.sub.c. Due to slotted resampling, the slots can be empty if no packet was transmitted sufficiently close to it.

[0097] The approach consists in setting a minimum number of sensing units per window, termed, that allows a sufficiently accurate reconstruction of the D signatures.

[0098] Then additional units are transmitted whenever the number of reflections of communication packets in the window is not sufficient to meet this minimum requirement.

[0099] The proposed method only requires the knowledge of whether a reflected communication packet is received in the current slot, i.e., no information about the future traffic pattern is needed.

TABLE-US-00002 Algorithm 2 - Injection of sensing units in window m. Input: M.sub.s. 1: # P1 - observation phase 2: N(m) no. of sensing units received in the first half of the window (either from reflected communications packets or injected). 3: # P2 - scheduling phase 4: N.sub.w (m) max(M.sub.s N.sub.a (m), 0). 5: Schedule S.sub.m = {s.sub.1, . . . , s.sub.Nw (m)}. 6: # P3 - transmission phase 7: for q =mW/2, . . . , (m + 1)W/2 1 do 8: if q S.sub.m then 9: if no reflected comm. packet received then 10: Transmit the sensing unit. 11: S.sub.m S.sub.m \ {q}. 12: else 13: Use the sensing unit from the comm. packet 14: S.sub.m S.sub.m \ {q}. 15: end if 16: else 17: if reflected comm. packet received then 18: Use the sensing unit from the comm. packet 19: S.sub.m S.sub.m \ {min.sub.sS }. 20: end if 21: end if 22: end for

[0100] The algorithm, summarized in the above Algorithm 2, operates in three phases, namely observation (P1), scheduling (P2) and transmission (P3). Recall that the D extraction already described follows a window-based approach, with subsequent windows overlapping by half of their length, as shown in FIG. 3.

[0101] Consider a time instant between the end of window m1 and the start of window m+1. This coincides with the half of window m, which is between slots mW/21 and mW/2. In this time instant it can be observed how many reflected communication packets were received in the first half of window m, which spans the indices from (m1)W/2 to mW/21 (P1, line 2 in Algorithm 2). This number is denoted as N.sub.0 (m). The whole injection algorithm operates on a half-window basis, as this allows reasoning on the sole current window m, as window m1 has ended and window m+1 has not yet started. Based on N.sub.Q (m), can computed how many sensing units would need in the remaining half of window m in order to meet the requirement of at least M.sub.s units, which is denoted by N.sub.w(m)=max(M.sub.sN.sub.a (m), 0). However, the sensing process has no knowledge of when future communication packets will be received, so the best that it can be done is to schedule the transmission of (m) sensing units in the next half-window. The slots in which these packets are scheduled can be selected according to a deterministic rule or a probability distribution. It is called S.sub.m={s.sub.1, . . . , s.sub.Nw(m)} the set of indices of the slots in which the additional sensing units for the next half-window (P2, lines 4-5 in Algorithm 2) are scheduled. While P1 and P2 are performed in a single time slot, before the second half of window m starts, P3 (lines 7-23 of Algorithm 2) is a dynamic process that spans the whole second half of window m. The indices of the slots considered in this part of the algorithm are q=mW/2, . . . , (m+1)W/21. Note that some unknown communication packets may be received in this second half-window. The procedure iterates over the slots and in each of them checks if a sensing unit was scheduled for that slot, i.e., if qS.sub.m. There are four possible cases: [0102] (1) qS.sub.m and no communication packet was received in this slot.

[0103] In this case the sensing unit is transmitted, q removed from Sm. [0104] (2) qS.sub.m and a communication packet (or more) was received in this slot. In this case the sensing unit is reused in the communication packet and q removed from Sm. [0105] (3) q.Math.S.sub.m and no communication packet was received in this slot.

[0106] In this case just move to the next slot without taking action. [0107] (4) q.Math.S.sub.m and a communication packet (or more) was received in this slot. In this case the sensing unit is reused in the communication packet, then the next sensing unit is removed from the scheduled ones, i.e., S.sub.m is set: S.sub.mS.sub.m \{mins.sub.sSm s}.

[0108] Note that, despite operating on a half-window basis, due to the overlap of adjacent windows, the algorithm only poses a constraint on the minimum number of packets sent per full window. This means that half a window can be empty as long as enough sensing units are received in the other half.

Scheduling the Sensing Units.

[0109] While P2, the scheduling of the sensing units, can be done with any arbitrary policy that guarantees that exactly N.sub.w(m) packets are scheduled in the next half window, the aim is to maximize the number of sensing units that can be piggybacked on communication packets, rather than using a dedicated transmission. From P3 in Algorithm 2, one can see that scheduling the sensing units towards the end of the half-window leaves more time for possible communication packets to become available and thus be reused instead of injecting a new sensing unit. Consequently, according to the present invention, the sensing units for the second half of window m as a burst of packets spaced by T, are scheduled, which occupy the last (m) slots of the window.

A Possible Exemplary Embodiment of the Computer Implemented Method

[0110] In the following, a possible embodiment of the computer implemented method of the invention is proposed, in view of the above description.

Step 1.

[0111] Select a slot duration T according to the desired maximum Doppler velocity of interest using the formula T=c/4 f.sub.ov.sub.max where v.sub.max is the desired maximum micro-Doppler velocity resolution in the spectrum, c is the light speed and f.sub.o is the carrier frequency.

Step 2.

[0112] Construct a sliding time window of length W slots. At each time frame, shift the window forward by an arbitrary frame duration, e.g., half of the window length (W/2).

Step 3.

[0113] The transmitter is used to transmit communication packets according to an underlying and ongoing communication process. The communication packets may contain CIR estimation fields, used by the processor to estimate the CIR. The communication can be performed using any communication protocol that supports, natively or via super resolution algorithms, a decimeter-level resolution in estimating the multiple paths in the CIR. These include mmWave communication systems such as IEEE 802.11ad/ay WLANs, 5G-NR (that naturally achieve high multipath estimation resolution) and sub 6 GHz communication standards such as IEEE 802.11a/g/n/ac (that can use super-resolution to achieve high multipath resolution). S4. Set an arbitrary integer threshold number of CIR samples needed in each window to obtain a good reconstruction of the CIR. The threshold can be set empirically or from theoretical bounds from sparse reconstruction theory.

Step 5.

[0114] At a pre-defined moment inside the window, if less than the threshold number of CIR estimates could be obtained from communication packets only, the processor schedules additional CIR estimation packets (containing only the CIR estimation field) to be transmitted in the current window, following an arbitrary policy, e.g., to transmit the K needed fields in the last K slots of the current window.

Step 6.

[0115] During the period of the current window that comes after the scheduling moment, the transmitter monitors how many communication packets containing CIR estimation fields are transmitted. For each such packet, the first additional remaining scheduled CIR estimation field is removed from the schedule as it is no longer needed to meet the threshold number of CIR measurements.

Step 7.

[0116] The transmitter transmits the scheduled CIR estimation fields. The receiver receives these additional waveforms in addition to the communication packets and the processor uses them to obtain additional CIR samples.

Step 8.

[0117] The processor the CIR estimates, h(t), computed from CIR estimation fields received in the current window (either from communication packets or from additionally injected CIR estimation fields). Denote the number of such estimates as N.sub.s, and each estimate is obtained at a sampling instant t.sub.i with i=1, . . . , N.sub.s. The CIR estimates are time-varying complex vectors of dimension L, where L is the number of different propagation paths in the signal.

Step 9.

[0118] The processor resamples the N.sub.s CIR estimates collected in the current window using the slotted resampling method. Starting from the available sequence h(t.sub.i) with i=1, . . . , N.sub.s, operate as follows: [0119] 1. Consider the current window as a regular grid with step size T and W total steps. [0120] 2. Construct a set of bins, each identified by index k, defined as the intervals [kTT/2, kT+T/2) with center kT. [0121] 3. Obtain the resampled sequence h(kT) using the following rule. For each bin k: [0122] a. If no sampling instant t.sub.i falls inside bin k, h(kT) is considered as a missing sample; [0123] b. If at least one instant t.sub.i falls inside bin k, set h(kT)=h(t.sub.k) where t.sub.k is the instant closest to the bin center kT among all sampling instants falling inside bin k. [0124] 4. Denote by U the number of available (non-missing) samples in the current window.

Step 10.

[0125] The processor constructs an U by W partial inverse Fourier matrix, F, as follows: [0126] 1. Build the W by W complete inverse Fourier matrix, whose element nm is equal to

[00016]$W^{-1/2}\exp \left(\frac{j2nm}{W}\right),$

where j is the imaginary unit. [0127] 2. Select, from the complete matrix, a U by W submatrix selecting the rows corresponding to the indices of the non-missing samples in the current window.

Step 11.

[0128] Select a subset of the paths where a target of interest is located, denoted by I1, I2, . . . , IQ identified by index q. This subset can be obtained using any side information such as a target tracking method.

Step 12.

[0129] Consider now independently each propagation path, identified by index by I1, I2, . . . , IQ. Denote by hq the U-dimensional complex vector of samples of the q-th selected path in the current window. Set up the following optimization problem

[00017]$\arg {\min }_{H_q}{.Math.h_q-{\mathrm{FH}}_q.Math.}_2^2\mathrm{such}\mathrm{that}{.Math.H_q.Math.}_0<S$

where H.sub.q is the Fourier transform of the (unknown) complete CIR path q in the current window (of dimension W), .sub.x is the x-norm of the argument vector and S is a parameter to be selected for the specific application.

Step 13.

[0130] Solve the optimization problem using any sparse optimization method, e.g., the iterative hard thresholding algorithm. Repeat S6-S7 for all paths I1, I2, . . . , IQ.

Step 14.

[0131] Compute the elementwise square magnitude (spectrum) of H.sub.q for the selected paths of interest.

Step 15.

[0132] Sum the spectra for the selected paths of interest, obtaining a total spectrum, H, for the current window.

Step 16.

[0133] Move forward the time window and repeat steps Step 3 to Step 10.

Step 17.

[0134] Stack together the spectrum vectors obtained from M adjacent windows to obtain the W by M micro-Doppler spectrogram.

A Possible Preferred Embodiment Implementation

[0135] In general terms, a system for joint communication and reconstruction of the micro-Doppler time-frequency spectrum from sparse channel measurements according to the invention comprises: [0136] a) a transmitter configured for transmitting wireless communication signals through a multipath channel, including channel estimation fields; [0137] b) a receiver, configured for receiving a wireless signal which is the reflection or refraction of the transmitted signal; and [0138] c) a processor, configured for implementing a method according to anyone of the possible embodiments described above.

[0139] Preferably, the transmitter and the receiver can be realized as a single transceiver sharing a same antenna array working in full-duplex mode.

[0140] According to possible embodiments, the estimated CIR is from a backscatter channel.

[0141] The transmitter can be equipped with an antenna array and phase shifters for directional beamforming.

[0142] Further, the CIR can be estimated for each of the different beampatterns used during the transmission.

[0143] In the following, it will be described one of the possible embodiments of the present invention implemented on a mmWave SDR platform. The implementation is based on the IEEE 802.11ay WiFi protocol, as it operates in the unlicensed 60 GHz band and supports CIR estimation for different BPs.

Testbed

[0144] The open-source mm-FLEX experimentation platform has been reproduced for a baseline design. It uses an FPGA-based baseband processor which can generate, capture and process (custom or standard compliant) frames with up to 1.76 GHz of bandwidth.

[0145] The baseband processor is connected to mmWave RF frontends with phased antenna arrays and supports various front-ends to operate in different frequency bands, e.g., at 28 GHz or 60 GHz. In the remainder of this description it will be used a 60 GHz RF front-end which simplifies experimentation as this is an unlicensed band, but it is noted that simply by changing the RF front-end, the present invention can operate in a different band, e.g., for 5G-NR compatibility.

[0146] In order to implement a system according to the present invention, the functionalities of the testbed have been augmented to enable full-duplex operation, not only in the baseband processor but also in the RF front-end. A block diagram of the system is shown in FIG. 4. In the baseband processor, the operation of the system is controlled by a state machine (SM) which triggers an AXI-DMA that reads the I/Q samples from the DDR-TX memory and feeds them to the DACs. The SM also triggers another AXI-DMA in the receiver datapath that saves the receives samples in the DDRRX memory, i.e., both datapaths are synchronized and no packet detection is required. To support the variable IFS extracted from real (or artificially generated) traces, a block RAM memory (BRAM) has been included in the FPGA logic that stores the IFS that will be used in the experiments. The SM reads these values sequentially, introducing a delay in the system according to the value read from memory. The variable IFS functionality can be disabled at runtime to configure a fixed IFS. It is remarked that since the up/down conversion stages are simultaneously used from the same mmWave development kit, the Tx/Rx sub-systems are fed by the same local oscillator and thus the Carrier Frequency Offset (CFO) is very low (<100 Hz), which enables the extraction of the D values required by the present invention.

[0147] IEEE 802.11ay CIR estimation details. In IEEE 802.11ay, in-packet beam tracking is introduced, where the CIR is estimated using different BPs within a single packet. This is done by appending a given number of training (TRN) fields to the packet. A TRN field is composed of 6 TRN units formed by complementary Golay sequences of 128 BPSK modulated samples, for a total of 768 samples.

TABLE-US-00003 TABLE 1 Summary of the implementation parameters. The suggested values based on experimental results are shown in bold. System parameters Grid step Tc 0.27 ms Window length W 64 Window overlap 32 Sparsity parameter {1, 2, 3, 4, 5, 6, 7} No. aggregated paths Q {1, 3, 5, 7, 9, 11, 13, 15} Min. no. measurements Ms {4, 8, 16, 24, 32, 64} IHT learning rate 1 IHT convergence threshold 10.sup.4 IHT maximum iteration number n.sub.max 200

[0148] According to the invention n.sub.TRN TRN fields are used as the sensing unit, where each TRN field employs a different BP, and n.sub.TRN is the number of subjects tracked by the system, as it suffices to use 1 TRN field per subject. Considering the typical number of people that are simultaneously tracked in human sensing systems, reasonable n.sub.TRN values go from 1 to 10. The CIR estimates obtained from each TRN field (BP) are then used as the inputs.

[0149] In Table 1 the system parameters used in the implementation are summarized. It is set T.sub.c=0.27 ms and W=64, which lead to (i) a velocity resolution of

[00018]$v=\frac{c}{2f_o{\mathrm{WT}}_c}0.14\mathrm{ms}$

and (ii) aliasing-free velocity measurements up to of

[00019]$v_{\mathrm{ma}x}=\frac{c}{4f_o{T}_c}4.48\text{m/s}.$

These values are not critical to the functioning of the system, and can be modified according to specific implementation requirements. However, for reliable D extraction without aliasing, it is advisable to adjust Tc to a value that allows capturing the range of velocities typically covered by human movement, e.g., approximately (2 to 3) m/s for a walking person, and up to 5 m/s for running or other fast movements. Note that suitable values of Tc can also be obtained in 5G-NR systems, where a base station can transmit downlink CSI-RS frames with a periodicity between 0.3125 ms and 80 ms. For a 5G-NR carrier frequency of 28 GHz, using T.sub.e=0.3125 ms leads to v.sub.max8.57 m/s, which is enough to capture fast human movement.

[0150] For people tracking, periodically transmitted in-packet beam training frames are used with twelve TRN units and antenna beams covering a FOV range from 45 to 45. Then, the already described distance and AoA estimation procedure are applied. Different values of M.sub.s, and Q are experimented, as reported in Table 1 and already described, while for the IHT algorithm the parameters that led to the most accurate convergence results on the experiments are selected, i.e., =1, =10-4 and n.sub.max=200.

[0151] In the following, the experimental results obtained with the above embodiment will be presented. The experiments were performed in a laboratory of 67 meters with a complex multi-path environment due to additional reflections caused by furniture, computers, screens, and a wide whiteboard.

Results on Synthetic Traces

[0152] As a first qualitative result, the D spectrograms obtained by the present invention on randomly sampled CIRs of a walking subject are shown in FIG. 5. For this, synthetic traces are used, generated by measuring the CIR using a uniform sampling interval equal to T.sub.c, and then setting to 0 a variable number of uniformly distributed values per window to simulate missing samples. This is a simplified case, as (i) the available (not removed) packets lie on a regular grid with spacing T.sub.c, therefore no approximation error is introduced by slotted resampling, and (ii) samples are removed on a per window basis, so a minimum number of packets in each window is guaranteed. Still, this evaluation is useful to highlight the impact of increasing the sparsity level of the measurements for the present invention compared to standard STFT. In the results here presented, no packet injection is performed, as the aim is to assess the impact of the number of measurements per window on the reconstructed D. In FIG. 5b the baseline walking spectrogram obtained using the standard STFT using the full window of 64 samples is shown.

[0153] The spectrogram shows a typical walking D modulation, with the contribution of the static clutter (the strong component at 0 velocity), of the torso (the strong oscillating component around 1.5 m/s and the limbs (the faint contributions around the torso component).

[0154] Moreover, a certain amount of noise and interference is present, as shown by the non-zero background level and the horizontal lines at around 2 m/s and 3.7 m/s. In FIG. 5c, the same method is applied to windows with only 16 out of 64 the samples retained, while the rest is set to 0. The impact is very strong as it completely corrupts the useful structure in the D signature. From FIG. 5d to FIG. 5f the results obtained by the present invention are shown, on the same sequence, with 64, 16 and 4 samples out of 64, respectively. Above each figure, it is also report the Root Mean-Squared Error (RMSE) of the D with respect to a ground truth spectrogram, shown in FIG. 5a. This was obtained from the STFT output with full measurement windows (FIG. 5b), by manually segmenting the useful D spectrum containing the gait information and setting to 0 all the background noise and the interference lines. Two interesting aspects are observed. On the one hand, the algorithm of the present invention can successfully recover the D spectrum even when a large fraction of the samples is missing, and the quality of the result decreases gradually with the sparsity of the available measurements. On the other hand, the present invention almost completely eliminates noise and interference that are present with standard STFT, thanks to the effect of the sparsity constraint in Eq. (7), obtaining a lower RMSE than STFT even when for the latter all measurements in the window are available. This aspect is the main reason why the present invention not only reduces the overhead needed for human sensing, but also improves its accuracy.

Realistic Traces: The Pdx/Vwave Dataset

[0155] Next, the performance of the present invention on realistic WiFi AP traces have been evaluated. This poses an experimental challenge, because commercial devices implementing the IEEE 802.11ay standard are not yet available, and no public datasets containing real traffic traces for the PHY layer of mmWave WiFi (IEEE 802.11ay/ad) exist. For this reason, the pdx/vwave dataset is used, containing real traffic traces captured in different real environments from WiFi APs employing a legacy (sub 6 GHz) WiFi protocol. Specifically, three, over one hour long, traces from this dataset are used, called psu cs, library and ug, respectively. Traces collected in different environments are selected to represent different kinds of traffic patterns (see Table 2 below).

TABLE-US-00004 TABLE 2 Details of the 3 sequences of the pdx/vwave dataset. Trace Environment No. frames Duration psu cs University CS dept. 260326 1:00 h library Public library 1300671 4:00 h ug Coffee shop 895721 2:34 h

[0156] The pdx/vwave dataset includes information about the transmission instants and packet sizes of all packets outgoing from the considered AP. Exploiting this information, measurements are performed transmitting packets according to these time patterns, using the BRAM in the FPGA to store the desired transmission instants. On top of the existing pdx/vwave communication patterns the injection algorithm (Algorithm 2) to send additional sensing units when needed is used.

[0157] Despite the fact that in the pdx/vwave dataset a legacy WiFi protocol is used, entailing much smaller packets sizes than IEEE 802.11ay, it is still reasonable to use it to obtain realistic packet transmission patterns. While in the pdx/vwave dataset the maximum physical layer PDU size is PPDUpdx=1.5 kB (without packet aggregation), in IEEE 802.11ay three main transmission modes are defined, namely High Throughput (HT), Directional Multi Gigabit (DMG) and Very High Throughput (VHT), with maximum physical layer PDU sizes, PPDUay, of 65 kB, 262 kB and 4692 kB, respectively. With the increase in the packet sizes, the data rates of mmWave systems have increased accordingly, and in IEEE 802.11ay they will range from 0.3 Gbps to several Gbps. As a numerical example, the traffic patterns in pdx/vwave with a typical bitrate of 4 Mbps would correspond to a bitrate of 0.7 Gbps in DMG IEEE 802.11 ay when using a packet size of 262 kB instead of 1.5 kB. Note that traces with a larger number of packets and smaller PDU sizes (as will likely be the case in real deployments) will simply increase the sensing accuracy and further reduce the overhead.

Human Activity Recognition Results

[0158] To evaluate the quality of the D spectrograms extracted by the present invention, such spectrograms are used as the input to a HAR method. Specifically, a standard approach is followed, training a deep neural network on a dataset of W dimensional D spectrograms, with =200 (equivalent to 1.76 s), in order to classify the movement performed by the person during that time. In order to provide a comparison with other IEEE 802.11ay HAR methods based on regular CIR sampling, such as RAPID, the following four activities are considered: walking, running, sitting and waving hands. For HAR, a standard Convolutional Neural Network (CNN) architecture is used, composed of four inception modules performing 11, 33 and 55 convolutions.

[0159] The number of filters used is 8, 16, 32 and 64 for the four modules, respectively. The convolutional blocks are followed by a fully-connected layer with 64 neurons, to which Dropout is applied, and a final Softmax layer with four outputs. The exponential-linear unit activation function is used after each layer.

Training Data

[0160] A training dataset has been collected involving six different subjects performing the four activities, for a total duration of about twelve minutes each. This leads to over 400 partially overlapping, 1.76 s long, D sequences per activity, which has been augmented as described shortly. Note that the training data only includes uniformly sampled CIR traces with sampling period Tc=0.27 ms. The CNN training is done for 80 epochs, using a learning rate of 10-4, the Adam optimizer and the cross-entropy loss function. In order to enhance the robustness of the CNN, an ad-hoc data augmentation strategy is applied: some of the CIR samples in each window of the training dataset are randomly removed, and then the IHT algorithm according to the invention is applied to reconstruct the spectrograms. The process is repeated using a sparsity level of , and , enlarging the training dataset to four times its original size, for a total of approximately 1600 D spectrograms per activity. A randomly selected subset of the training data (around 10%) was used as a validation set to tune the hyperparameters of the CNN.

Test Data

[0161] The CNN has been tested on the D spectrograms obtained from CIR samples collected using the pdx/vwave packet traces previously described. Four, randomly selected, 20 s long traces (one per activity) are collected for each of the three sequences types (psu cs, library and ug). The experiments are repeated for different values of the minimum number of sensing units per window, M.sub.s=4, 8, 16, 32, 64, for a total of sixty test sequences. The test data involves a single subject, which was not included in the training set.

HAR F1 Score.

[0162] The performance of the CNN is evaluated with the per-class F1 score metric, which effectively summarizes the precision and recall and preserves the class-specific results. FIG. 7 shows the total average per-class F1 score over the 60 sequences, for different values of the minimum number of sensing units per window, M.sub.s. As a baseline for comparison, the F1 score obtained by the RAPID algorithm it is also reported, which extracts the D signatures by regularly sampling the CIR. The results show that the present invention can reach over 0.9 F1 scores on all activities with M.sub.s=8 already, which corresponds to only of the full measurements window. Notably, with M.sub.s=4, the low number of measurements per window affects significantly only the Sitting and Waving hands activities, which involve fine-grained movements and are therefore more difficult to classify. Finally, the results from the present invention and RAPID are compared. For a fair comparison, RAPID's STFT has been implemented to extract the D and trained the CNN on the resulting spectrograms without enlarging the dataset using different levels of sparsity described in the previous section. Instead, the training procedure of is directly used, since it has been found that the sparsity-based data augmentation slightly reduced RAPID's performance. It can be seen that sparse recovery problem formulation and enforcing a sparsity constraint on the individual paths is beneficial to HAR performance. The gap is particularly significant for Sitting and Waving hands as they involve lower energy traces in the spectrograms; these are more easily corrupted by noise and interference, that the present invention is mostly able to reject (see, again, the comparison between FIG. 5b and FIG. 5d).

Overhead Analysis

[0163] Increasing M.sub.s to improve the HAR performance also increases the overhead of the present invention. A first general measure of this can be obtained comparing the maximum size of a PPDU in IEEE 802.11 ay to the size of a sensing unit. Recalling the three different modes previously described and the size of an IEEE 802.11ay TRN field (768 bits), it is obtained that a sensing unit, with nTRN=1, is 0.1%, 0.03% and 0.002% of a PPDU in HT, DMG and VHT, respectively. Moreover, the channel occupation time for a sensing unit with nTRN=1 is 436 ns, which is a negligible fraction (0.16%) of a slot of duration T.sub.c.

[0164] Next, to evaluate the overhead of the present invention on a realistic communication scenario, the traces of the pdx/vwave dataset are used. In this way, it can also be assessed the impact of injecting sensing units, as they are not useful to the communication process. Denote by ci the number of bits in the i-th communication packet transmitted in the trace. As the number of bits transmitted in each trace refers to a legacy, lower bitrate, WiFi protocol, the packet sizes according to the maximum PHY layer packet size in IEEE 802.11 ay are rescaled. The size of packet i in each trace as

[00020]${\tilde{c}}_i=\left(\frac{\mathrm{PPDUay}}{\mathrm{PPDUpdx}}\right)*c_i$

are rescaled, with PPDUay=262 kB. It is called n.sub.c the number of transmitted communication packets in a trace, by TRN.sub.len the length, in bits, of a piggybacked or injected TRN field, by n.sub.inj the number of injected sensing units and by nTRN the number of TRN fields used in every sensing operation (it is considered fixed, whereas in reality it is determined by the number of subjects in the environment). The overhead is defined as a function of n.sub.inj as

[00021]$\begin{array}{cc}\mathrm{OH}\left(n_{\mathrm{inj}}\right)=\frac{n_{\mathrm{TRN}}\left(n_C+n_{\mathrm{inj}}\right){\mathrm{TRN}}_{\mathrm{len}}}{{.Math.}_{i=1}^{n_c}{\tilde{c}}_i}& \left(10\right)\end{array}$

[0165] In FIG. 8A the overhead obtained on each of the three pdx/vwave traces is shown, using nTRN=1. The overhead for different values of nTRN can be obtained by using it as a multiplicative factor on the values in FIG. 8A. It is noticed that the overhead scales almost linearly as M.sub.s is increased from 4 to 64. For values of M.sub.s<32, the entailed overhead is less than 4%, falling below 1% for M.sub.s=8 As a reference, the overhead for M.sub.s=64 is reported, which is the value obtained by injecting sensing units continuously into the channel, piggybacking them eventually on communication packets if possible. Note that existing approaches requiring uniform CIR sampling, like RAPID, would require an even higher overhead, as not only do they need 64 samples per window, but these samples have to be regularly spaced as no resampling procedure is carried out. This means they would have to take precedence over potential data packets so that they are sent exactly at the right sampling time.

[0166] From FIG. 8B, one can see that the present invention can achieve an F1 score of over 0.9 for every activity for a minimum of M.sub.s=8 sensing units per window, resulting in a sensing overhead of less than 1%. With this configuration, the present invention achieves a better F1 score than existing approaches, while reducing overhead by a factor of seven and being compatible with random access MAC protocols.

Sensitivity to the Choice of the Parameters

[0167] In FIG. 9, the effect of varying parameters Q it is shown, representing the number of paths aggregated around the person's position (see Eq. (9)), and , which is the maximum number of resolvable Doppler components, equal to the sparsity parameter in the IHT algorithm.

[0168] The HAR per-class F1 score are computed using a random subset of the 60 test sequences. The values adopted in the experiments are reported in Table 1.

Impact of Changing Q.

[0169] The results show that the present invention is robust to almost any value of Q when considering walking and running, whereas sitting and waving hands are negatively affected by reducing Q below seven. This is due to the fact that while walking and running are, in most cases, distinguishable even from the sole contribution of the torso, this is not true for sitting and waving hands that require including the reflection paths coming from the limbs. Computational complexity considerations are also in order for high values of Q, as it leads to solving Q times Eq. (7) at each D extraction process. As the problems are independent, they can be solved in parallel, and thus a reasonable approach is to tune Q according to a trade-off between D reconstruction accuracy and hardware resource availability for parallelization. In the following, Q=9 is used. Considering that B=1.76 GHz transmission bandwidth (one IEEE 802.11 ay channel) is used, the range resolution is c/2B=8.5 cm. This means that summing the contribution of Q/2 paths before and after the one corresponding to the torso, a region of 34 cm around the person's position is included in the spectrum, which is a reasonable value considering typical body sizes and that the subjects are moving.

[0170] Fixing Q=9, in FIG. 9B, it is shown that the best values for are 2 and 3 for all the activities.

[0171] This is because using =1 often leads to only reconstructing the 0 Doppler component in the spectrogram, losing the information on the person's movement. On the other hand, choosing (too high makes the IHT reconstruction imprecise, as with a low number of measurements per window enforcing more sparsity is beneficial to restrict the number of possible solutions to Eq. (7).

[0172] In this description, the first ISAC system that can sense human D signatures from irregular and sparse CIR estimates has been designed and implemented. These estimates are obtained in a standard compliant way by both reusing optional CIR estimation fields appended to communication packets and by the sporadic injection of sensing packets whenever communication traffic is absent. Different from the existing ISAC methods, the present invention is based on a sparse recovery approach to the D reconstruction, which is theoretically grounded in the instrinsic sparse multi-path environment of the channel. This enables an accurate D extraction from a significantly lower number of randomly distributed CIR samples, thus drastically reducing the sensing overhead. After a CIR resampling step along the time domain, the present invention performs an iterative sparse reconstruction in the frequency domain, decoupling different propagation paths at first, to leverage their sparsity property, and then combining them to obtain the final D spectrum.

[0173] While the present invention is compatible with different communication systems (e.g., 3GPP 5G-NR, and IEEE WLANs), for the implementation an IEEE 802.11ay SDR platform working in the 60 GHz band has been used. The system has been tested on a large set of standard-compliant CIR traces matching the traffic patterns of real WiFi access points, performing a typical downstream application such as HAR. The results show that the present invention entails over seven times lower overhead compared to prior methods, while achieving better performance.

[0174] The present invention has been so far described with reference to preferred embodiments thereof. It is to be meant that each one of the technical solutions implemented in the preferred embodiments, herein described by way of example, can advantageously be combined, differently from what described, with the other ones, to create additional embodiments, belonging to the same inventive core and however all within the protective scope of the here below reported claims.